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Abstract 

We review the current status of studies of the coalescence of binary neutron star systems. 
We begin with a discussion of the formation channels of merging binaries and we discuss the 
most recent theoretical predictions for merger rates. Next, we turn to the quasi-equilibrium 
formalisms that are used to study binaries prior to the merger phase and to generate initial 
data for fully dynamical simulations. The quasi-equilibrium approximation has played a key 
role in developing our understanding of the physics of binary coalescence and, in particular, 
of the orbital instability processes that can drive binaries to merger at the end of their life- 
times. We then turn to the numerical techniques used in dynamical simulations, including 
relativistic formalisms, (magneto-)hydrodynamics, gravitational-wave extraction techniques, 
and nuclear microphysics treatments. This is followed by a summary of the simulations per- 
formed across the field to date, including the most recent results from both fully relativistic 
and microphysically detailed simulations. Finally, we discuss the likely directions for the field 
as we transition from the first to the second generation of gravitational-wave interferometers 
and while supercomputers reach the petascale frontier. 
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1 Introduction 



Binaries composed of neutron stars (NSs) and black holes (BHs) have long been of interest to astro- 
physicists. They provide many important constraints for models of massive star evolution and com- 
pact object formation, and are among the leading potential sources for detection by gravitational- 
wave (GW) observatories. While it remains uncertain whether mergers of compact binaries are 
an important contributor to the production of r-process elements, they are now thought to be the 
leading candidate to explain short-duration, hard-spectrum gamma-ray bursts (often abbreviated 
to "short-hard" GRBs, or merely SGRBs). 

The first neutron star-neutron star (NS-NS) binary to be observed was PSR B1913-f 16, in 
which a radio pulsar was found to be in close orbit around another NS |134| . In the decades 
since its discovery, the decay of the orbit of PSR B1913-I-16 at exactly the rate predicted by 
Einstein's general theory of relativity (see, e.g., [3011 1320] ) has provided strong indirect evidence 
that gravitational radiation exists and is indeed correctly described by general relativity (GR). 
This measurement led to the 1993 Nobel Prize in physics for Hulse and Taylor. 

According to the lowest-order dissipative contribution from GR, which arises at the 2.5PN level 
(post-Newtonian; where the digit indicates the expansion order in [w/c]^ in the Taylor expansion 
term), and assuming that both NSs may be approximated as point masses, a circular binary orbit 
decays at a rate da/dt = —a/TQw where a is the binary separation and the gravitational radiation 
merger timescale tqw is given by 
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where Mi, M2, and M = AIi + M2 are the individual NS masses and the total mass of the 
binary, /i = M1M2/M is the reduced mass, q = M2/M1 is the binary mass ratio, and we assume 
geometrized units where G = c = 1 (as we do throughout this paper, unless otherwise noted). The 
timescale for an elliptical orbit is shorter, and it can be shown that eccentricity is reduced over 
time by GW emission, leading to a circularization of orbits as they decay. A quick integration 
shows that the time until merger is given by Tmcrgo — tgw/4. 
The luminosity of such systems in gravitational radiation is 
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which, at the end of a binary's lifetime, when the components have approached to within a few 
NS radii of each other, is comparable to the luminosity of all the visible matter in the universe 
(~ 10^'^ erg/s). The resulting strain amplitude observed at a distance D from the source (assumed 
to be oriented face-on) is given approximately by 
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The first measurement that wiU hkely be made with direct GW observations is the orbital decay 
rate, with the period evolving (for the circular case) according to the relation 
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where T is the orbital period and uj the angular frequency, and thus the "chirp mass. 
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is likely to be the easiest parameter to determine from GW observations. 

Several NS-NS systems are now known, including PSR J0737-3039 [SS], a binary consisting of 
two observed pulsars, which allows for the prospect of even more stringent tests of GR |148j . Even 
with the handful of observed sources to date, one may use this sample to place empirical limits on 
the expected rate of NS-NS mergers [142] and to constrain the many parameters that enter into 
population synthesis calculations |213j . With regard to the former, the very short merger timescale 
for J0737, Tmergo =85 Myr, makes it especially important for estimating the overall rate of NS-NS 
mergers since it is a priori very unlikely to detect a system with such a short lifetime. 

Although black hole-neutron star (BH-NS) binaries are expected to form through the same 
processes as NS-NS binaries, none has been detected to date. This is generally thought to reflect 
their lower probability of detection in current surveys, in addition to intrinsically smaller numbers 
compared to NS-NS systems [39]. BH-NS systems are an expected byproduct of binary stellar 
evolution, and properties of the population may be inferred from population synthesis studies 
calibrated to the observed NS-NS sample (see, e.g., [57]). 

In this review, we will summarize the current state of research on relativistic mergers, beginning 
in Section [2] with a description of the astrophysical processes that produce merging binaries and 
the expected parameters of these systems. The phases of the merger are briefly described in 
Section [3] In Section [31 we discuss the numerical techniques used to generate quasi-equilibrium 
(QE) sequences of NS-NS configurations, and we summarize the QE calculations that have been 
performed. These sequences yield a lot of information about NS physics, particularly with regard 
to the nuclear matter equation of state (EOS). They also serve as initial data for dynamical merger 
calculations, which we discuss next, focusing in turn on the numerical hydrodynamics techniques 
used to compute mergers and the large body of results that has been generated, in Sections [S] and 
[SI respectively. We pay particular attention to how numerical studies have taken steps toward 
answering a number of questions about the expected GW and electromagnetic (EM) emission from 
merging binaries, and we discuss briefly the possibility that they may be the progenitors of SGRBs 
and a source of r-process elements. We close with conclusions and a look to the future in Section [71 

While most of this review focuses on NS-NS mergers, many of the methods used to study NS-NS 
binaries are also used to evolve BH-NS binaries, and it has become clear that both merger types 
may produce similar observational signatures as well. For a review focusing on BH-NS merger 
calculations, we encourage the reader to consult the recent work by Shibata and Taniguchi 279^ . 
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2 Evolutionary Channels and Population Estimates 



Merging NS-NS and BH-NS binaries, i.e., those for which the merger timescale is smaller than the 
Hubble time, are typically formed through similar evolutionary channels in stellar field populations 
of galaxies |37] (both may also be formed through dynamical processes in the high-density cores of 
some star clusters, but the overall populations are smaller and more poorly constrained; see |254| 
for a review). It is difficult to describe the evolutionary pathways that form NS-NS binaries 
without discussing BH-NS binaries as well, and it is important to note that the joint distribution 
of parameters such as merger rates and component masses that we could derive from simultaneous 
GW and EM observations will constrain the underlying physics of binary stellar evolution much 
more tightly than observing either source alone. 

Population synthesis calculations for both merging NS-NS and BH-NS binaries typically favor 
the standard channel in which the first-born compact object goes through a common-envelope 
(CE) phase, although other models have been proposed, including recent ones where the progenitor 
binary is assumed to have very nearly-equal mass components that leave the main sequence and 
enter a CE phase prior to either undergoing a supernova [42 , 54 . Simulations of this latter process 
have shown that close NS-NS systems could indeed be produced by twin giant stars with core 
masses > 0.15 Mq, though twin main sequence stars typically merge during the contact phase [174] . 

In the standard channel (see, e.g., [44111771. and Figure [1] for an illustration of the process), the 
progenitor system is a high- mass binary (with both stars of mass M > 8 — 10 Mq to ensure a pair 
of supernovae). The more massive primary evolves over just a few million years before it leaves 
the main sequence, passes through its giant phase, and undergoes a Type lb, Ic, or II supernova, 
leaving behind what will become the heavier compact object (CO): the BH in a BH-NS binary 
or the more massive NS in an NS-NS binary. The secondary then evolves off the main sequence 
in turn, triggering a CE phase when it reaches the giant phase and overflows its Roche lobe. 
Dynamical friction shrinks the binary separation dramatically, until sufficient energy is released 
to expel the envelope. Without this step, binaries would remain too wide to merge through the 
emission of GWs within a Hubble time. Eventually, the exposed. Helium-rich core of the secondary 
undergoes a supernova, either unbinding the system or leaving behind a tight binary, depending 
on the magnitude and orientation of the supernova kick. 

This evolutionary pathway has important effects on the physical parameters of NS-NS and 
BH-NS binaries, leading to preferred regions in phase space. The primary, which can accrete some 
matter during the CE phase, or during an episode of stable mass transfer from the companion 
Helium star, should be spun up to rapid rotation (see |177| for a review). In NS-NS binaries, we 
expect that this process will also reduce the magnetic field of the primary down to levels seen in 
"recycled" pulsars, typically up to four orders of magnitude lower than for young pulsars [1791 [75] . 
The secondary NS, which never undergoes accretion, is likely to spin down rapidly from its nascent 
value, but is likelier to maintain a stronger magnetic field. 

While this evolutionary scenario has been well studied for several decades, many aspects remain 
highly uncertain. In particular; 

• The CE efficiency, which helps to determine the expected range of binary separations and 
the mass of the primary compact object after its accretion phase, remains very poorly con- 
strained ( 2251 [571 [751 1333] . If too much matter is accreted by the NS, it may undergo 
accretion- induced collapse to a BH though the growing body of observed NS-NS sys- 
tems docs help place constraints on the allowed range of accretion- related parameters. 

• The exact relation between a star's initial mass and the eventual compact object mass is bet- 
ter understood, but significant theoretical uncertainties remain, and the relation is sensitive 
to the metallicity (often unknown), mainly through the effects of mass loss in stellar winds, 
and to the details of the explosion itself [3311 1201] . 
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Figure 1: Cartoon showing standard formation channels for close NS-NS binaries through binary 
stellar evolution, taken from |177J . 

• The maximum allowed NS mass will affect whether the primary remains a NS or undergoes 
accretion-induced collapse to a BH; its value is dependent upon the as-yet undetermined 
nuclear matter EOS. At present, the strongest limit is set by the binary millisecond pulsar 
PSR J1614-2230, for which a mass of A/ns = 1.97± 0.04 Mq was determined by Shapiro time 
delay measurements |81l . Determining the NS EOS from GW observations may eventually 
provide stronger constraints }233| I217[ I182j , including a determination of whether supernova 
remnants are indeed classical hadronic NS or instead have cores consisting of some form of 
strange quark matter or other elementary particle condensate [2191 11181 12261 11191 [231 IS] • 



The supernova kick velocity distribution is only partially understood, leading to uncertainties 
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as to which systems wih become unbound after the second explosion. [132| I319| 12151 1151) 



Given all these uncertainties, it is reassuring that most estimates of the NS-NS and BH-NS 
merger rate, expressed either as a rate of mergers per Myr per "Milky Way equivalent galaxy" or 
as a predicted detection rate for LIGO (the Laser Interferometer Gravitational- Wave Observatory) 
and Virgo (see Section ISTSl below) . agree to within 1-2 orders of magnitude, which is comparable to 
the typical uncertainties that remain once all possible sources of error are folded into a population 
synthesis model. In Table [1] we show the predicted detection rates of NS-NS and BH-NS mergers 
for both the first generation LIGO detectors ( "LIGO" ) , which ran at essentially their design spec- 
ifications [2], and the Advanced LIGO ("AdLIGO") configuration due to go onhne in 2015 [287] . 
We note that the methods used to generate these results varied widely. In [142j . the authors used 
the observed parameters of close binary pulsar systems to estimate the Galactic NS-NS merger 
rate empirically (such results do not constrain the BH-NS merger rate) . In jl94l 1126) , the two 
groups independently estimated the binary merger rate from the observed statistics of SGRBs. 
In these cases, one does not get an independent prediction for the NS-NS and BH-NS merger 
rate, but rather some linear combination of the two. In both cases, the authors estimated that, if 
NS-NS and BH-NS mergers are roughly equal contributors to the observed SGRB sample, LIGO 
will detect about an order of magnitude more BH-NS mergers since their higher mass allows them 
to be seen over a much larger volume of the Universe. As they both noted, should either type 
of system dominate the SGRB sample, we would expect a doubling of LIGO detections for that 
class, and lose our ability to constrain the rate of the other using this method. Many population 
synthesis models have attempted to understand binary evolution within our Galaxy by starting 
from a basic parameter survey of the various assumptions made about CE evolution, supernova 
kick distributions, and other free parameters. In [3181179) . population synthesis models are normal- 
ized by estimates of the star formation history of the Milky Way. In |139l 1214) , parameter choices 
are judged based on their ability to reproduce the observed Galactic binary pulsar sample, which 
allows posterior probabilities to be applied to each model in a Bayesian framework. A review by 
the LIGO collaboration of this issue may be found in 1 . 

Should the next generation of GW interferometers begin to detect a statistically significant 
number of merger events including NSs, it should be possible to constrain several astrophysical 
parameters describing binary evolution much more accurately. These include 

• The relative numbers of BH-NS and NS-NS mergers: Interferometric detections are sensitive 
to a binary's "chirp mass" (see Eq. |6]), and to the binary mass ratio as well [56 | ITT | l3T6| 1315) 
if the signal-to-noise ratio is sufficiently high. Even if the merger signal takes place at 
frequencies too high to fall within the LIGO band, it should still be possible in most cases 
to determine whether the primary's mass exceeds the maximum mass of a NS. 

• The mass ratio probability distribution for merging binary systems: If both binary compo- 
nents' masses are determined, we will be able to constrain both the BH mass distribution 
in merging binaries and the NS mass ratio distribution. Knowledge of the former would 
determine, e.g., whether the current low estimates for the mass accreted onto the primary 
CO core during the CE phase are correct [38], as this model predicts that BH masses in 
close BH-NS binaries should cluster narrowly around A/bh = W Mq. Previous calculations 
assuming larger accreted masses typically favored mass ratios closer to unity, since the pri- 
maries often began as NSs and underwent "accretion-induced collapse" to a BH during the 
CE phase. The NS-NS mass distribution will allow for tests of whether "twins", i.e., systems 
whose zero-age main sequence (or "ZAMS") masses are so close that they both leave the 
main sequence before either undergoes a supernova, play a significant role in the merging 
NS-NS population 05]. 

Interestingly enough, while it has generally been assumed that Advanced LIGO or another 
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Table 1: Estimated initial and advanced LIGO rates for BH-NS and NS-NS mergers from pop- 
ulation synthesis calculations and other methods. The methods used are, in order, empirical 
constraints from the observed sample of binary pulsars ('Empirical'), constraints on the combined 
NS-NS/BH-NS merger rate assuming that they are the progenitors of short- hard gamma-ray bursts 
('SGRBs'), population synthesis models calibrated to the star formation rate in the Milky Way 
('Pop. Synth. - SFR'), and population synthesis calibrated against the observed Galactic binary 
pulsar sample ('Pop. Synth. - NS-NS'). We note that observations of binary pulsars do not yield 
constraints for BH-NS binaries. SGRB observations may produce constraints on NS-NS merger 
rates, BH-NS, merger rates, or both, depending on which sources are the true progenitors, but 
this remains unclear. Therefore, the table quotes results assuming a roughly equal split between 
the two. The official review of these results and their implications by the LIGO/Virgo Scientific 
Collaborations may be found in [T]. 
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interferometer will make the first direct observations of NS-NS mergers and their immediate after- 
math, such systems may have already been observed and simply not interpreted correctly. Indeed, 
while virtually all NS-NS merger calculations performed to date that consider EM emission have 
looked at the high-energy prompt emission, Nakar and Piran suggest that the mass ejection from 
mergers should yield an observable radio afterglow |196) . Moreover, they argue that such post- 
merger afterglows may have already been observed as a ~month-long radio burst occurring in the 
outer regions of a galaxy. While such an outburst could also result from a supernova, the lumi- 
nosity required would be an order of magnitude larger than those previously observed. Given the 
length and timescales characterizing radio bursts, no NS-NS simulation has been able to address 
the model directly, but it certainly seems plausible that the time-variable magnetic fields within 
a stable hypermassive remnant could generate enough EM energy to power the resulting radio 
burst [275) . These potential EM observations of mergers are likely to spur further research into the 
amount and velocity of merger ejecta, which could then be coupled to a larger-scale astrophysical 
simulation of the potential radio afterglow. 
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3 Stages of a Binary Merger 



The qualitative evolution of NS-NS mergers, or indeed any compact binary merger, has long been 
understood, and may be divided roughly into inspiral, merger, and ringdown phases, each of which 
presents a distinct set of challenges for numerical modeling and detection. As a visual aid, we 
include a cartoon summary in Figure [2l originally intended to describe black hole - black hole 
(BH-BH) mergers, and attributed to Kip Thorne. Drawn before the advent of the supercomputer 
simulations it envisions, we note that merger waveforms for all compact binary mergers have proven 
to be much smoother and simpler than shown here. To adapt it to NS-NS mergers, we note that 
NSs are generally assumed to be essentially nonspinning, and that the "ringdown" phase may 
describe either a newly formed BH or a NS that survives against gravitational collapse. 




*on owrt ^ a u perooTY-tputet^ known' 



Figure 2: Cartoon picture of a compact binary coalescence, drawn for a BH-BH merger but 
applicable to NS-NS mergers as well (although NSs are generally assumed to be nonspinning). 

Summarizing the evolution of the system through the three phases: 

• After a pair of supernovae yields a relatively tight NS-NS binary, the orbital separation decays 
over long timescales through GW emission, a phase that takes up virtually all of the lifetime 
of the binary except the last few milliseconds. During the inspiral phase, binary systems 
may be accurately described by QE formalisms, up until the point where the gravitational 
radiation timescale becomes comparable to the dynamical timescale. The evolution in time 
is well-described by FN expansions, currently including all terms to 3FN [47 , though small 
deviations can arise because of finite-size effects, especially at smaller separations (see Eqs.[T] 
-|6]for the lowest-order 2.5FN expressions for circular inspirals). 

• Once the binary separation becomes no more than a few times the radii of the two NSs, bi- 
naries rapidly become unstable. The stars plunge together, following the onset of dynamical 
instability, and enter the merger phase, requiring full GR simulations to understand the com- 
plicated hydrodynamics that ensues. According to all simulations to date, if the NS masses 
are nearly equal, the merger resembles a slow collision, while if the primary is substantially 
more massive than the secondary the latter will be tidally disrupted during the plunge and 
will essentially accrete onto the primary. Since the NSs are most likely irrotational just prior 
to merging, there is a substantial velocity discontinuity at the surface of contact, leading to 
rapid production of vortices. Meanwhile, some fraction of the mass may be lost through the 
outer Lagrange points of the system to form a disk around the central remnant. This phase 
yields the maximum GW amplitude predicted by numerical simulations, but with a signal 
much simpler and more quasi-periodic than in the original cartoon version. GW emission 
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during the merger encodes important information about the NS EOS, particularly the GW 
frequency /cut at which the binary orbit becomes unstable (see Eq. H]) resulting in a character- 
istic cutoff in GW emission at those frequencies. Meanwhile, the merger itself may generate 
the thermal energy that eventually powers a SGRB, which occurs when the neutrinos and 
anti-neutrinos produced by shock-heated material annihilate around the remnant to produce 
high-energy photons. 

• Finally, the system will eventually settle into a new, dynamically stable configuration through 
a phase of ringdown, with a particular form for the GW signal that depends on the remnant's 
mass and rotational profile. If the remnant is massive enough, it will be gravitationally un- 
stable and collapse promptly to form a spinning BH. Otherwise it must fall into one of three 
classes depending on its total mass. Should the remnant mass be less than the maximum 
mass Miso supported by the nuclear matter EOS for an isolated, non-rotating configuration, 
it will clearly remain stable forever. Instead a remnant that is "supramassive," i.e., with a 
mass above the isolated stationary mass limit but below that allowed for a uniformly rotating 
NS (typically < 1.2Miso, with very weak dependence on the EOS; see, e.g., |146l[7n[72] . and 
references therein) may become unstable. Supramassive remnants are stable against gravita- 
tional collapse unless angular momentum losses, either via pulsar- type emission or magnetic 
coupling to the outer disk, can drive the angular velocity below the critical value for stability. 
If the remnant has a mass above the supramassive limit, it may fall into the hypermassive 
regime, where it is supported against gravitational collapse by rapid differential rotation. 
Hypermassive NS (HMNS) remnants can have significantly larger masses, depending on the 
EOS (see, e.g., [30 [ I270| I85 | I288| I113j ). and will survive for timescales much longer than the 
dynamical time, undergoing a wide variety of oscillation modes. They can emit GWs should 
a triaxial configuration yield a significant quadrupole moment, and potentially eject mass 
into a disk around the remnant. Eventually, some combination of radiation reaction and 
magnetic and viscous dissipation will dampen the differential rotation and lead the HMNS to 
collapse to a spinning BH, again with the possibility that it may be surrounded by a massive 
disk that could eventually accrete. The energy release during HMNS collapse provides the 
possibility for a "delayed" SGRB, in which the peak GW emission occurs during the initial 
merger event, but the gamma-ray emission, powered by the collapse of the HMNS to a BH, 
occurs significantly later. 

Most calculations indicate that a geometrically thick, lower-density, gravitationally bound 
disk of material will surround whatever remnant is formed. Such disks, which are geomet- 
rically thick, are widely referred to as "tori" throughout the literature, though there is no 
clear distinction between the two terms, and we will use "disk" throughout this paper to 
describe generically the bound material outside a central merger remnant. Such disks are 
expected to heat up significantly, and much of the material will eventually accrete onto the 
central remnant, possibly yielding observable EM emission. Given the low densities and rel- 
atively axisymmetric configuration expected, disks are not significant GW emitters. There 
may be gravitationally unbound outfiow from mergers as well, though dynamical simulations 
neither confirm nor deny this possibility yet. Such outfiows, which can be the sites of exotic 
nuclear reactions, are frequently discussed in the context of r-process element formation, but 
their inherently low densities make them difficult phenomena to model numerically with high 
accuracy. 



3.1 Comparison to BH— NS mergers 

Since this three-phase picture is applicable to BH-NS mergers as well, it is worthwhile to compare 
the two merger processes at a qualitative level to understand the key similarities and differences. 
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Inspiral for BH-NS mergers is also well-described by PN expansions up until shortly before the 
merger, but the parameter space is fundamentally different. First, since BHs are heavier than 
NSs, the dynamics can be quite different. Also, since BHs may be rapidly-spinning (i.e., have 
dimensionless spin angular momenta as large as J/M^ 1), spin-orbit couplings can play a very 
important role in the orbital dynamics of the binary, imprinting a large number of oscillation modes 
into the GW signal (see, e.g., |124[ [56]). From a practical standpoint, the onset of instability in 
BH-NS mergers should be easier to detect for LIGO, Virgo, and other interferometers, since the 
larger mass of BH-NS binaries implies that instability occurs at lower GW frequencies (see Eq. 21 
noting that the separation a at which mass-transfer begins scales roughly proportionally to the BH 
mass). 

The onset of instability of a BH-NS binary is determined by the interplay of the binary mass- 
ratio, NS compactness, and BH spin, with the first of these playing the largest role (see Figures IS- 
IS of |296| and the summary in .279J). In general, systems with high BH masses and/or more 
compact NSs tend to reach a minimum in the binding energy as the radius increases, leading 
to a dynamical orbital instability that occurs near the classical innermost stable circular orbit 
(ISCO). In these cases, the NS plunges toward the BH before the onset of tidal disruption, and is 
typically "swallowed whole" . leaving behind almost no material to form a disk. The GW emission 
from such systems is sharply curtailed after the merger event, yielding only a low- amplitude, 
high-frequency, rapidly-decaying "ringdown" signal (see, e.g., [153] ). Numerical calculations have 
shown that even in borderline cases between dynamical instability and mass-shedding the NS is 
essentially swallowed whole, especially in cases where the BH in either non-spinning or spinning in 
the retrograde direction, which pushes the ISCO out to larger radii (see, e.g., [285 ,.284|I278[ [9T 1 [93 ] ). 

A richer set of phenomena occurs when the BH-NS mass ratio is closer to unity, the NS is 
less compact, the BH has a prograde spin direction, or more generally, some combination of those 
factors. In such cases, the NS will reach the mass-shedding limit prior to dynamical instability, and 
be tidally disrupted. Unlike what was described in semi-analytic Newtonian models (see, e.g., [69l 
I224| [77] ) and seen in some early Newtonian and quasi- relativistic simulations (see, e.g., |1641 
I165[ 1137] . stable mass transfer, in which angular momentum transfer via mass-shedding halts the 
inspiral, has never been seen in full GR calculations, nor even in approximate GR models (see 
the discussion in W). Even so, unstable mass transfer can produce a substantial disk around the 
BH, though in every GR simulation to data the substantial majority of the NS matter is accreted 
promptly by the BH (see |279) for a detailed summary of all current results). The exact structure 
of the disk and its projected lifetime depend on the binary system parameters, with the binary 
mass ratio and spin both important in determining the disk mass and the BH spin orientation 
critical for determining both the disk's vertical structure and its thermodynamic state given the 
shock heating that occurs during the NS disruption. In general, the mass and temperature of 
the post-merger disks are comparable to those seen in some NS-NS mergers, and inasmuch as 
either is a plausible SGRB progenitor candidate then both need to be viewed as such. To date, no 
calculation performed in full GR has found any bound and self-gravitating NS remnant left over 
after the merger, including both NS cores that survive the initial onset of mass transfer by recoiling 
outward (predicted for cases in which stable mass transfer was thought possible, as noted above) 
or those in which bound objects form through fragmentation of the circum-BH disk. Motivated 
by observations of wide-ranging timescales for X-ray fiares in both long and short GRBs 221], 
the latter channel has been suggested to occur in collapsars [223] and mentioned in the context 
of mergers [239] . possibly on longer timescales than current simulations permit. Even so, there is 
no analogue to the HMNS state that may result from NS-NS mergers, nor any mechanism for a 
delayed SGRB as provided by HMNS collapse. 
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3.2 Qualitative numerical results 

Constructing QE sequences for a given set of NS parameters requires sophisticated numerical 
schemes, but not supercomputer-scale resources, as we discuss in Section |4] below, focusing first 
on the numerical techniques used to construct QE binary data in GR, and the astrophysical 
information contained in the GW emission during the inspiral phase. Merger and ringdown, on 
the other hand, typically require large-scale numerical simulations, including some of the largest 
calculations performed at major supercomputer centers, as we discuss in detail Section [S] and [B] 
below. 

To illustrate the various physical processes that occur during NS-NS mergers, we show the 
evolution of three different NS-NS merger simulations in FiguresOIH and[5l taken from Figures 4- 
6 of |143j . In Figure 131 we see the merger of two equal-mass NSs, each of mass Mns = 1.4 M©, 
described by the APR (Akmal-Pandharipande-Ravenhall) EOS [3]. In the second panel, clear 
evidence of "tidal lags" is visible shortly after first contact, leading to a slightly off-center collision 
pattern. By the third panel, an ellipsoidal HMNS has been formed, surrounded by a disk of 
material of lower density, which gradually relaxes to form a more equilibrated HMNS by the final 
panel. In Figure SI we see a merger of two slightly heavier equal-mass NSs with M-^g = 1.5 Mq. 
In this case, the deeper gravitational potential limits the amount of mass that goes into the disk, 
and once a BH is formed (with a horizon indicated by the dashed blue circle in the final panel) 
it accretes virtually all of the rest mass initially present in the two NSs, with only 0.004% of the 
total remaining outside the horizon. 

In Figure m we see the merger of an unequal- mass binary, with masses Mi = 1.3 Mq and 
M2 = 1.6 Af0. In this case, the heavier NS partially disrupts the lighter NS prior to merger, 
leading to the secondary NS being accreted onto the primary. In this case, a much more massive 
disk is formed and, even after a BH forms in the center of the remnant, a substantial amount of 
matter, representing 0.85% of the total mass, remains outside the horizon. Later accretion of this 
material could potentially release the energy required to power a SGRB. 
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Figure 3: Isodensity contours and velocity profile in the equatorial plane for a merger of two 
equal-mass NSs with Mns = 1.4 M© assumed to follow the APR model [3] for the NS EOS, taken 
from Figure 4 of [143] . The hypermassive merger remnant survives until the end of the numerical 
simulation. 

4 Initial Data and Quasi-Equilibrium Results 
4.1 Overview 

While dynamical calculations are required to understand the GW and EM emission from BH- 
NS and NS-NS mergers, some of the main qualitative features of the signals may be derived 
directly from QE sequences. From the variation of total system energy with binary angular velocity 
along a given sequence, it is possible to construct an approximate GW energy spectrum dEQ^/df 
immediately from QE results, essentially by performing a numerical derivative (see FigurelHl). Doing 
so for a number of different sequences makes it possible to identify key frequencies where tidal effects 
may become measurable and to identify these with binary parameters such as the system mass 
ratio and NS radius. Similarly, since QE sequences should indicate whether a binary begins to shed 
mass prior to passage through the ISCO (see Figure[7]), one may be able to classify observed signals 
into mass-shedding and non-mass-shedding events, and to use the critical point dividing those cases 
to help constrain the NS EOS. Single-parameter estimates have been derived for NS-NS binaries 
using QE sequences [98] (and for BH-NS binaries using QE j295] and dynamical calculations [278j ). 
NS-NS binaries typically approach instability at frequencies few ^ 1 kHz, where laser shot noise 
is severely degrading the sensitivity of an interferometer detector. To observe ISCO-related effects 
with higher signal-to-noise, it may be necessary to operate GW observatories using narrow-band 
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Figure 4: Isodensity contours and velocity profile in the equatorial plane for a merger of two equal- 
mass NSs with Mns = 1.5 M© assumed to follow the APR model [3^ for the NS EOS, taken from 
Figure 5 of [143) . With a higher mass than the remnant shown in Figure [31 the remnant depicted 
here collapses promptly to form a BH, its horizon shown by the dashed blue circle, absorbing all 
but 0.004% of the total rest mass from the original system. 



signal recycling modes, in which the sensitivity in a narrow range of frequencies is enhanced at the 
cost of lower sensitivity to broadband signals [ST] . 

It is important to note that, while the potential parameter space for NS EOS models is still 
very large, a much smaller set may serve to classify models for comparison with the first generation 
of GW detections. Indeed, by breaking up the EOS into piecewise polytropic segments, one may 
use as few as four parameters to roughly approximate all known EOS models, including standard 
nuclear models as well as models with kaon or other condensates |233) . To illustrate this, we show 
in Figure[5]four different QE models for NS-NS configurations with different EOS, taken from j300] : 
all have Mi = 1.15 Mq and M2 = 1.55 Mq, and they correspond to the closest separation for which 
the QE code still finds a convergent result. 

The inspiral of NS-NS binaries may yield complementary information about the NS structure 
beyond what can be gleaned from QE studies of tidal disruption. NSs have a wide variety of 
oscillation modes, including f-modes, g-modes, and r-modes, any of which may be excited by 
resonances with the orbital frequency as the latter sweeps upward. Should a particular oscillation 
mode be excited resonantly, it can then serve briefiy as an energy sink for the system, potentially 
changing the phase evolution of the binary. For example, in a rapidly spinning NS, excitation of 
the m — I r-mode can be significant, yielding a change of over 100 radians for the pre-merger GW 
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Figure 5: Isodensity contours and velocity profile in the equatorial plane for a merger of two 
unequal-mass NSs with Mi = 1.3 M© and M2 = 1.6 M©, with both assumed to follow the APR 
model [3 for the NS EOS, taken from Figure 6 of |143j . In unequal-mass mergers, the lower mass 
NS is tidally disrupted during the merger, forming a disk-like structure around the heavier NS. 
In this case, the total mass of the remnant is sufhciently high for prompt collapse to a BH, but 
0.85% of the total mass remains outside the BH horizon at the end of the simulation, which is 
substantially larger than for equal-mass mergers with prompt collapse (see Figure |4]). 



signal phase in the case of a millisecond spin period |160) . For NS-NS mergers in the field, this 
would require one of the NSs to be a young pulsar that has not yet spun down significantly, which 
is unlikely because of the difficulty in obtaining such an extremely small binary separation after 
the second supernova. Other modes, such as the / = 2 f-mode, may be excited in less extreme 
circumstances, also yielding information about NS structure parameters |105j . 

4.2 Quasi-equilibrium formalisms 

It has long been known that the GW emission from eccentric binaries is very efficient at radiating 
away angular momentum relative to the radiated energy |222j : as a result, the orbital eccentricity 
decreases as a binary inspirals, so that orbits should be very nearly circular long before they 
enter the detection range of ground-based interferometers. The only exception could be from a 
dynamical capture process that would create a binary with a significant eccentricity and very small 
orbital separation. Such eccentric binaries have been predicted to form in the nuclear cluster of 
our Galaxy (see, e.g., j209) ) and in core-collapsed globular clusters |125lll66j . However, at present, 
no formalism exists to construct initial data for such systems, besides superposing the individual 
components with sufficiently large initial separations to minimize constraint violations. 

Using this circularity of primordial binaries as a starting point, one may use the constraint 
equations of GR, along with an assumption of quasi-circularity, to derive sets of elliptic equations 
describing compact binary configurations. For both QE and dynamical calculations, most groups 
typically make use of the Arnowitt-Deser-Misner (ADM) 3+1 splitting of the metric which 
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Figure 6: Diniensionless binding energy Eb/Mo vs. dimensionless orbital frequency A/oil, where 
Mo is the total ADM (Arnowitt-Deser-Misner) mass of the two components at infinite separation, 
for two QE NS-NS sequences that assume a piecewise polytropic NS EOS, taken from Figure 16 
of [300] . The equal- mass case assumes Mns = 1-35 Mq for both NSs, while the unequal- mass case 
assumes Mi = 1.15 and M2 = 1.55 Mq. The thick curves are the numerical results, while the 
thin curves show the results from the 3PN approximation. The lack of any minimum suggests that 
instability for these configurations occurs at the onset of mass shedding, and not through a secular 
orbital instability. 



foliates the metric into a set of three-dimensional hypersurfaces by introducing a time coordinate. 
The resulting form of the metric, which is completely general, is written 

gf,„ = (-a^ + PiP'-^dt^ + 2/3jdt dx' + jijdx'dx', (7) 

where a is known as the lapse function, /3i the shift vector, and the spatial three-metric intrinsic 
to the hypersurface. We are following the standard relativistic notation here where Greek indices 
correspond to four-dimensional quantities and Latin indices to spatial three-dimensional quantities. 
Thus, the shift vector is a 3-vector, raised and lowered with the spatial 3-metric 7^ rather than 
the spacetime 4-metric g^^. To simplify matters, one typically defines a conformal factor -0 that 
factors out the determinant of the 3-metric, such that 

^^[det(7,,)]'^''- (8) 

introducing the conformal 3-metric jij = ip^^^ij with unit determinant. While the 3-metric is a 
fundamental component of the geometric structure of the spacetime, the lapse function and shift 
vector are gauge quantities that simply refiect our choice of coordinates. Thus, while one often 
determines the lapse and shift in order to construct a appropriately "stationary" solution in the 
relevant coordinates between neighboring time slices, their values are often replaced to initialize 
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Figure 7: Mass-shedding indicator x = f^^g^ ) / ( ^^5^ ) orbital frequency MqQ, where 

V ^ / cq V ^ / pole 

h is the fluid enthalpy and the derivative is measured at the NS surface in the equatorial plane 
toward the companion and toward the pole in the direction of the angular momentum vector, for 
a series of QE NS-NS sequences assuming equal-mass components, taken from Figure 19 of [300] . 
Here x = 1 corresponds to a spherical NS, while x = indicates the onset of mass shedding. More 
massive NSs are more compact, and thus able to reach smaller separations and higher angular 
frequencies before mass shedding gets underway. 



dynamical runs with more convenient choices and thus different assumptions about coordinate 
evolution in time. 

The field equations of general relativity take the deceptively simple form 

Gfii, = Rf_ii, — -g^ii^R — SttT^i, (9) 

where G^i, is the Einstein tensor, R^,^ and R the Ricci curvature tensor and the curvature scalar, 
and Tfin the stress-energy tensor that accounts for the presence of matter, electromagnetic fields, 
and other physical effects that contribute to the mass-energy of the spacetime. Since GR is a 
second-order formulation, valid initial data must include not only the metric but also its first time 
derivative. It generally proves most convenient to introduce the time derivative of the metric after 
subtracting away the Lie derivative with respect to the shift, yielding a quantity known as the 
extrinsic curvature, Kij: 

{dt-Cp)^,,=-2aK,,. (10) 

Both the 3-metric and extrinsic curvature are symmetric tensors with six free parameters. 

For systems containing NSs, one must consider the effects of nuclear matter through its presence 
in the stress-energy tensor T'^'^. It is common to assume that the matter has the EOS describing 
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Figure 8: Isodensity contours for QE models of NS-NS binaries, taken from Figures 9-12 of [300] . 
In each case, the two NSs have masses Mi = 1.15 Mq (left) and M2 = 1.55 Mq (right), and 
the center-of-mass separation is as small as the QE numerical methods allow while able to find a 
convergent result. The models assume different EOS, resulting in different central concentrations 
and tidal deformations. 



a perfect, isotropic fluid, for which the stress energy tensor is given by 

T^'' = (po + Po£ + n^i^^i'' + (11) 

where po, e, P and u'^ are the fluid's rest-mass density, speciflc internal energy, pressure, and 
4-velocity, respectively. Many calculations further assume that the NS EOS is described by an 
adiabatic polytrope, for which 

P=(r-l)poe = fcpr, (12) 

where F is the adiabatic index of the gas and k a constant, though a number of models designed 
to incorporate nuclear physics and/or strange matter condensates have also been constructed and 
studied (see Sections 14.41 and [6] below) . 

The problem in constructing initial data is not so much producing solutions that are self- 
consistent within GR, but rather to specify a sufficient number of assumptions to fully constrain 
a solution. Indeed, there are only four constraints imposed by the equations of GR, known as 
the Hamiltonian and momentum constraints. The Hamiltonian constraint is found by projecting 
Einstein's equations twice along the direction defined by a normal observer, and describes the way 
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stress-energy leads to curvature in the metric (see, e.g., [28] for a thorough review): 

R + K^ - KijK'^ = 167rp, (13) 

where R is the scalar curvature of the 4-mctric, K ~ Kl is the trace of the extrinsic curvature, and 

p = n • T • n = a^yO" = ph{au°f - P (14) 

is the total energy density seen by a normal observer. The third term indicates that the total 
energy density is found by projecting the stress-energy tensor in the direction of the unit-length 
timelike normal vector n, whose components are given by 

= (-a, 0,0,0). (15) 

In the final expression h = 1 + e + P/ po is the specific enthalpy of the fluid, and the combination 
r„ = au'^ represents the Lorentz factor of the matter seen by an inertial observer. The notation 
here makes use of the standard summation convention, in which repeated indices are summed over. 

Projecting Einstein's equations in the space and time directions leads to the vectorial momen- 
tum constraint 

D,K} - D,K = 8^ji, (16) 

where Di represents a three-dimensional covariant derivative and ji = pohTnUi is the total mo- 
mentum seen by a normal observer. 



4.2.1 The Conformal Thin Sandwich formalism 

In order to specify all the free variables that remain once the Hamiltonian and momentum con- 
straints are satisfied, two different techniques have been widely employed throughout the numerical 
relativity community. One, known as the Conformal Transverse- Traceless (CTT) decomposition, 
underlies the Bowen-York [52] solution for black holes with known spin and/or linear momentum 
that is widely used in the "moving puncture" approach. To date, however, the CTT formalism has 
not been used to generate NS-NS initial data, and we refer readers to [2791 [M] for descriptions of 
the CTT formalism applied to BH-NS and BH-BH initial data, respectively. 

To date, most groups have used the Conformal Thin Sandwich (CTS) formalism to generate 
QE NS-NS data (see [55] for a review, [T3] 1136] for the initial steps in the formulation, and |321[ 
13221 13281 170] for derivations of the form in which it is typically used today) . One first specifies that 
the conformal 3-metric is spatially flat, i.e., 7y = Sij, where Sij is the Kronecker delta function. 
Under this assumption, the only remaining parameter defining the spatial metric is the conformal 
factor 'ijj, which serves the role of a gravitational potential. Indeed, in the limit of weak sources, 
it is linearly related to the standard Newtonian potential. Next, one specifies that there exists a 
helical Killing vector, so that, as the configuration advances forward in time, all quantities remain 
unchanged when properly rotated at constant angular velocity in the azimuthal direction. This is 
sufficient to fix all but the trace of the extrinsic curvature, with the other components forced to 
satisfy the relation 



2aV''^ 



(17) 



The trace of the extrinsic curvature K remains a free parameter in this approach. While one 
may choose arbitrary prescriptions to fix it, most implementations choose a maximal slicing of the 
spatial hypersurfaces by setting K — dtK = 0. Under these assumptions the Hamiltonian and 
momentum constraints, along with the trace of the Einstein equations, yield five elliptic equations 
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for the lapse, shift vector, and confornial factor: 

= -V'' ( ^K'^K.j - 27rp ) , (18) 



72 ' ' N ' ' 



(a^) = a^p'^i-ij^KijK'^ +2Tr^'^{p + 2S), (19) 



= 2aV'^if^JVj(aV'"^) + 16aV'V, (20) 



1, 
3 

where 

5 = (g^. + n^n,)r^'' = fi-j = 3P + (p + P) [l - T'^] (21) 
is the trace of the stress-energy tensor projected twice in the spatial direction. While these five 
equations are linked and the right-hand sides are non-linear, they are amenable to solution using 
iterative methods. Boundary conditions are set by assuming asymptotic fiatness: at large radii, 
the metric takes on the Minkowski form so a — > 1, — > 1, and /S^^^ — > x r . We note that a purely 
corotating shift term yields zero when we apply the left-hand side of Eq. [20l so we may subtract 
it away and solve the equation with a boundary condition of zero instead. 

The breakdown in Eqs.[T8l[T9l and [20] is not unique. The Meudon group |1231I176] . to pick one 
example, has often chosen to define = Ina and /3 = ln(a^'^), and replace Eqs. 1181 and \W\ with 
the equivalent pair 

W^i^ = ij^K'^K.j -Wii'W'l3 + 4TrtP'^(p + S), 
4 ^ 2 

This approach is sufficient to define the field component of the configuration, but one still 
needs to solve for the matter quantities as well. One starts by assuming that there is a known 
prescription for reconstructing the density, internal energy, and pressure from the enthalpy h. 
Next, one has to assume some model for the spin of the NS. While corotation is often a simpler 
choice, since the velocity field of the matter is zero in the corotating frame, the more physically 
reasonable condition is irrotational flow. Indeed a realistic NS viscosity is never suflaciently large 
to tidally lock the NS to its companion during inspiral [45] 1145] . If we define the co- momentum 
vector Wi — hui, irrotational flow implies the vanishing of its curl: 

V X w — V^Wjy — Vi/W^ — d^w^ — d^w^ — 0, (22) 

which allows us to deflne a velocity potential ^ such that w = V^*. Using these quantities, one 
may write down the integrated Euler equation 

/iar„/(l -7yC/*C/^') = const., (23) 

where the 3-velocity C/* of the fluid with respect to an Eulerian observer is given by 

u'.^^. (24) 

and the orbital 3-velocity C/q with respect to the same observer by C/q = For details on the 

ways in which one may construct an elliptic equation for the velocity potential, we refer to the 
derivation in [123j . 

To date, all QE sequences and dynamical runs in the literature have assumed that NSs are 
either irrotational or synchronized, but it is possible to construct the equations for arbitrary NS 
spins so long as they are aligned [52] 1304] . While suggestions are also given there on how to 
construct QE sequences with intermediate spins using the new formalism, none have yet appeared 
in the literature. Similarly, a formalism to add magnetic fields self-consistently to QE sequences 
has been constructed [309] . as current dynamical simulations typically begin from data assuming 
either zero magnetic fields or those that only contribute via magnetic pressure. 
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4.2.2 Other formalisms 



The primary drawback of tlie CTS system is the lack of generality in assuming the spatial metric to 
be conformally flat, which introduces several problems. The Kerr metric, for example, is known not 
to be conformally flat, and conformally fiat attempts to model Kerr BHs inevitably include spurious 
GW content. The same problem affects binary initial data: in order to achieve a configuration 
that is instantaneously time-symmetric, one actually introduces spurious gravitational radiation 
into the system, which can affect both the measured parameters of the initial system as well as 
any resulting evolution. 

Other numerical formalisms to specify initial data configurations in GR have been derived using 
different assumptions. Usui and collaborators derived an elliptic set of equations by allowing the 
azimuthal component of the 3-metric to independently vary from the radial and longitudinal com- 
ponents |314[|313j . finding good agreement with the other methods discussed above. A number of 
techniques have been developed to construct helically symmetric spacetimes in which one actually 
solves Einstein's equations to evaluate the non-conformally-fiat component of the metric, which 
are typically referred to as "waveless" or "WL" formalisms [2561 [501 1286) . In terms of the funda- 
mental variables, rather than specifying the components of the conformal spatial metric by ansatz, 
one specifies instead the time derivative of the extrinsic curvature using a physically motivated 
prescription. These methods are designed to match the proper asymptotic behavior of the metric 
at large distances, and may be combined with techniques designed to enforce helical symmetry of 
the metric and gauge in the near zone (the near zone helical symmetry, or "NHS" formalism) to 
produce a global solution |310|I329|I311) . QE sequences generated using this formalism |310j have 
shown that the resulting conformal metric is indeed non-fiat, with deviations of approximately 1% 
for the metric components, and similar differences in the system's binding energy when compared 
to equivalent CTS results. They suggest |311| that underestimates in the quadrupole deformations 
of NS prior to merger may result in total phase accumulation errors of a full cycle, especially for 
more compact NS models. 

QE formalisms reflect the assumption that binaries will be very nearly circular, since GW emis- 
sion acting over very long tiniescales damps orbital eccentricity to negligible values for primordial 
NS-NS binaries between their formation and final merger. Binaries formed by tidal capture and 
other dynamical processes, which may be created with much smaller initial separations, are more 
likely to maintain significant eccentricities all the way to merger (see, e.g., 12091 for a discussion of 
such processes for BH-BH binaries) and it has been suggested based on simple analytical models 
that such mergers, likely occurring in or near dense star clusters, may account for a significant 
fraction of the observed SGRB sample 166 . However, more detailed modeling is required to work 
out accurate estimates of merger rates given the complex interplay between dynamics and binary 
star evolution that determines the evolution of dense star clusters, and given the large uncertainties 
in the distributions of star cluster properties in galaxies throughout the universe. No initial data 
have ever been constructed in full GR for merging NS-NS binaries with eccentric orbits since the 
systems are then highly time-dependent, while the calculations performed to evolve them generally 
use a superposition of two stationary NS configurations jl21) . 

4.3 Numerical implementations 

There are a number of numerical techniques that have been used to solve these elliptic systems. The 
first calculations of NS-NS QE sequences, in both cases for synchronized binaries, were performed 
by Wilson, Mathews, and Marronetti [3221 1323] [187] and Baumgarte et al. [25l|26]. The former used 
a finite differencing scheme, and centered different quantities at cells, vertices, and faces in order 
to construct a system of equations that was solved using fast matrix inversion techniques, while 
the latter used a Cartesian multigrid scheme, restricted to an octant to increase computational 
efficiency. 
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After a formalism for evaluating QE irrotational NS-NS sequences was developed [51] I308j . 
some of the first results were obtained by Uryii and Eriguchi, who developed a finite-differencing 
code in spherical coordinates allowing for the solution of relativistic NS-NS binaries using Green's 
functions |3081 1312] , Their method extended the self-consistent field (SCF) work of |146] . which 
had previously been applied to axisymmetric configurations. Irrotational configurations were also 
generated by Marronetti et al. |183j . using the same finite difference scheme as found in the work 
on synchronized binaries. 

The most widely used direct grid-based solver in numerical relativity is the Bam_Elliptic 
solver |55| . which solves elliptic equations on single rectangular grids or multigrid configurations. 
It is included within the Cactus code, which is widely used in 3-D numerical relativity |302J . In 
particular it has been used to initiate a number of single and binary BH simulations, including one 
of the original breakthrough binary puncture works [61 . 

Lagrangian methods, typically based on smoothed particle hydrodynamics (SPH) [1801 11171 
I190J ) have been used to generate both synchronized and irrotational configurations for PN [101 l99l 
fTUn [Too] and conformally flat (CF) [Ml Ml HQH [2031 HQS [2041 [34l ESI [33] calculations of 
NS-NS mergers, but they have not yet been extended to fully GR calculations, in part because of 
the difficulties in evolving the global spacetime metric. 

The most widely used data for numerical calculations are those generated by the Meudon 
group (see Section below for details on their calculations and [123) for a detailed description 
of their methods) . The code they developed, Lorene |176| , uses multidomain spectral methods to 
solve elliptic equations (while the code has been used primarily for relativistic stellar and binary 
configurations, it can be used as a more general solver). Around each star, one creates a set of 
nested, contiguous grids, with points arrayed in the radial and angular directions. The innermost 
grid has spheroidal geometry, and the surrounding grids are annular. The outermost grid may 
be allowed to extend to spatial infinity through a compactification transformation of the radial 
coordinate. To solve elliptic equations for various field quantities, one breaks each into a sum of 
two components, each of whose source terms are concentrated in one NS or the other. Similarly, 
the source terms themselves are split into two pieces, ideally so each component is well-described 
by spheroidal spectral coefficients centered around each star. Using the spectral expansion, one 
may pass values from one star to the other and then recalculate spectral coefficients for the other 
grid configuration. This scheme has several efficiency advantages over direct grid-based methods, 
which helps to explain its popularity. First, the domain geometry may be chosen to fit to a NS 
surface, which eliminates Gibbs phenomenon-related errors and allows for exponential convergence 
with respect to the number of grid points, rather than the geometric convergence that characterizes 
finite difference-based grid codes. Second, the use of spectral methods requires much less computer 
memory than grid-based codes, and, as a result, Lorene is a serial code that can run easily on any 
off-the-shelf PC, rather than requiring a supercomputer platform. 

4.4 Quasi-equilibrium and pre-merger simulations 

NS-NS binaries may be well approximated by QE configurations up until they reach separations 
comparable to the sizes of the binary components themselves, that latter phase lasting a fraction 
of a second after an inspiral of millions of years or more. The eventual merger will occur after the 
binary undergoes one of two possible orbital instabilities. If the total binary energy and angular 
momenta reach a minimum at some separation, which defines the ISCO, the binary becomes 
dynamically unstable and plunges toward merger. Alternately, if the NS fills its Roche Lobe 
(typically the lower density NS) mass will transfer onto the primary and the secondary will be 
tidally disrupted. The parameters of some NS-NS systems could technically allow for stable 
mass transfer, in which mass loss from a lighter object to a heavier one leads to a widening of the 
binary separation. This does occur for some binaries containing white dwarfs, but every dynamical 
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calculation to date using full GR or even approximate GR has found that the rapid inspiral rate 
leads to inevitably unstable mass transfer and the prompt merger of a binary. 

Many of the results later confirmed using relativistic QE sequences were originally derived in 
Newtonian and PN calculations, particularly as explicit extensions of Chandrasekhar's body of work 
(see [66]). Chandrasekhar's studies of incompressible fluids were first extended to compressible 
binaries by Lai, Rasio, and Shapiro [155| I154[ I157| I156| I158J . who used an energy variational 
method with an ellipsoidal treatment for polytropic NSs. They established, among other results, 
the magnitude of the rapid inspiral velocity near the dynamical stability limit |155j , the existence 
of a critical polytropic index (n « 2) separating binary sequences undergoing the two different 
terminal instabilities jl54j . the role played by the NS spin and viscosity and magnitude of finite- 
size effects in relation to IPN terms |157|I156] . and the development of tidal lag angles as the binary 
approaches merger |158J . They also determined that for most reasonable EOS models and non- 
extreme mass ratios, as would pertain to NS-NS mergers, an energy minimum is inevitably reached 
before the onset of mass transfer through Roche lobe overflow. The general results found in those 
works were later confirmed by '197', who used a SCF technique |129II130] . finding similar locations 
for instability points as a function of the adiabatic index of polytropes, but a small positive offset 
in the radius at which instability occurred. Similar results were also found by |306| I307j , but with 
a slight modification in the total system energy and decrease in the orbital frequency at the onset 
of instability. 

The first PN ellipsoidal treatments were developed by Shibata and collaborators using self- 
consistent fields [263 EMI Elil [276| [295 and by Lombardi, Rasio, and Shapiro fT75^. Both groups 
found that the nonlinear gravitational effects imply a decrease in the orbital separation (increase 
in the orbital frequency) at the instability point for more compact NS. This result reflects a fairly 
universal principle in relativistic binary simulations: as gravitational formalisms incorporate more 
relativistic effects, moving from Newtonian gravity to IPN and on to CF approximations and finally 
full GR, the strength of the gravitational interaction inevitably becomes stronger. The effects seen 
in fully dynamical calculations will be discussed in Section [6l below. 

The first fully relativistic GTS QE data for synchronized NS-NS binaries were constructed by 
Baumgarte et al. [211 ES], using a grid-based elliptic solver. Their results demonstrated that the 
maximum allowed mass of NSs in close binaries was larger than that of isolated NSs with the 
same (polytropic) EOS, clearly disfavoring the "star-crushing" scenario that had been suggested 
by [322| I185j using a similar GTS formalism (but see also the error in these latter works addressed 
in |104j . discussed in Section[63]below). Baumgarte et al. also identified how varying the NS radius 
affects the ISGO frequency, and thus might be constrained by GW observations. Using a multigrid 
method. Miller et al. jlSSj showed that while confornial fiatness remained valid until relatively near 
the ISGO, the assumption of sycnronized rotation broke down much earlier. Usui et al. |314) used 
the Green's function approach with a slightly different formalism to compute relativistic sequences 
and determined that the GTS conditions were valid up until extremely relativistic binaries were 
considered. 

The first relativistic models of physically realistic irrotational NS-NS binaries were constructed 
by the Meudon group [51] using the Lorene multi-domain pseudo-spectral method code. Since 
then, the Meudon group and collaborators have constructed a wide array of NS-NS initial data, 
including polytropic NS models ^123, ,298, ,299^ , as well as physically motivated NS EOS mod- 
els [36] or quark matter condensates |168) . Irrotational models have also been constructed by 
Uryii and collaborators [308[ 1312) for use in dynamical calculations, and nuclear/quark matter 
configurations have been generated by Oechslin and collaborators |208l 1205) . A large compila- 
tion of QE GTS sequences constructed using physically motivated EOS models including EPS 
(Friedman-Pandharipande) |218) . SLy (Skyrme Lyon) f83^, and APR [3] models, along with piece- 
wise polytropes designed to model more general potential cases (see [233]), was published in 300]. 

The most extensive set of results calculated using the waveless/near-zone helical symmetry 
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condition appear in |311l, with equal-mass NS-NS binary models constructed for the FPS, Sly, 
and APR EOS in addition to F = 3 polytropes. Results spanning all of these QE techniques are 
summarized in Table [51 
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Table 2: A summary of various studies focusing on QE sequences of NS-NS binaries. Please refer 
to Section [5] for a discussion of papers that focus on dynamical simulations instead. Gravita- 
tional schemes include Newtonian gravity ('Newt.'), lowest-order post-Newtonian theory ('PN'), 
confornial thin sandwich ('CTS') including modified forms of the spatial metric ('Mod. CTS'), 
and waveless/near-zone helical symmetry techniques. Numerical methods include ellipsoidal for- 
malisms ('EUips.'), self-consistent fields ('SCF'), numerical grids ('Grid'), multigrids, and multi- 
patch, Green's function techniques ('Green's'), spectral methods ('Spectral'), or SPH relaxation 
('SPH'). With regard to EOS models, 'WD' refers to the exact white dwarf EOS assuming a cold 
degenerate electron gas [65 . The 'Physical' EOS models include the EPS |218| . SLy [83], and 
APR [3] nuclear EOS models, along with their parameterized approximations and other physically 
motivated models. The compactness C = M/R refers to the value for a NS in isolation before it is 
placed in a binary, and plays no role in Newtonian physics. The mass ratio q — M2/M1 is defined 
to be less than unity, and 'spin' refers to either synchronized or irrotational configurations. 
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Lai 
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Newt. 


Ellips. 


r = |,|,2,oo 
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1.0 




Syn. 
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Newt. 
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Syn. 


Shibata 
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Shibata 
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PN 
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1.0 




Syn. 


Lombardi 


|175| 


PN 
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1.0 




Syn/Irr. 


Taniguchi 
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Limousin 
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CTS 


Spectral 


Quark 
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Oechslin 
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Irr. 
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5 Dynamical Calculations: Numerical Techniques 



5.1 Overview: General relativistic (magneto-) hydrodynamics and mi- 
crophysical treatments 

NS-NS binaries are highly relativistic systems, numerous groups now run codes that evolve both 
GR metric fields and fluids self-consistently, with some groups also incorporating an ideal magne- 
tohydrodynamic evolution scheme that assumes infinite conductivity. The codes that evolve the 
GR hydrodynamics or magnetohydrodynamics (GRHD and GRMHD, respectively) equations are 
many and varied, incorporating different spatial meshes, relativistic formalisms, and numerical 
techniques, and we will summarize the leading variants here. All full GR codes now make use of 
the significant insight gained from BH-BH merger calculations, but much work on these systems 
predates the three 2005 "breakthrough" papers by Pretorius [227) . Goddard [22], and the group 
then at UT Brownsville (now at RIT) [61], with the first successful NS-NS merger calculations 
announced already in 1999 [ 282j . A list of the groups that have performed NS-NS merger cal- 
culations using full GR is presented below; note that many of these groups have also performed 
BH-NS simulations, as discussed in the review by Shibata and Taniguchi |279j . 

Of the full GR codes used to evolve NS-NS binaries, almost all are grid-based and make use of 
some form of adaptive mesh refinement. The one exception is the SpEC code developed by the SXS 
collaboration, formed originally by Caltech and Cornell, which has used a hybrid spectral-method 
field solver with grid-based hydrodynamics. Most make use of the BSSN formalism for evolving 
Einstein's equations (see Sec. 15.2.11 below), while the HAD code uses the alternate Generalized 
Harmonic Gauge (GHG) approach. This technique is also used by the SXS collaboration and the 
Princeton group, who have both performed simulations of merging BH-NS binaries (see Sec. 16.61) 
but have yet to report any results on NS-NS mergers. Three groups have reported results for 
NS-NS mergers including MHD (HAD, Whisky, and UIUC), while the KT (Kyoto/Tokyo) group 
has reported magnetized evolutions of HMNS remnants (see |275] and references therein for a 
discussion of their work and that of other numerical relativity groups), but have yet to use that 
code for a NS-NS merger calculation. 

While full GR codes were being developed to study NS-NS binaries, a parallel and rather 
independent track developed to study detailed microphysical effects in binary mergers using ap- 
proximate relativistic schemes. This includes codes like that developed by the MPG group that 
accurately track the production of neutrinos and antineutrinos and their annihilation during a 
merger, as well as post-processing routines that use extensive nuclear chains to track the produc- 
tion of rare high-atomic number r-process elements in merger ejecta [122J . Meanwhile, the Bremen 
group's SPH code includes variable-temperature physically motivated equations of state |243| and 
magnetohydrodynamics [229], and has been used with a multi-group fiux-limited diffusion neutrino 
code to generate expected neutrino signatures from merger calculations [82] . A summary of groups 
performing NS-NS merger calculations is presented in Table |31 

5.2 GR Numerical techniques 
5.2.1 GR formalisms and gauge choice 

There are two distinct schemes used in all binary merger calculations performed to date, the BSSN 
(Baumgarte-Shapiro-Shibata-Nakamura) (2721 [27] and Generalized Harmonic formalisms. For 
general reviews of these formalisms, as well as other developments in numerical relativity, we refer 
readers to two recent books on numerical relativity [U [31] • Here we only briefly summarize these 
two schemes. 

The BSSN formalism was adapted from the 3+1 ADM approach, with quantities deflned as in 
Eqs. [7] and [U While the original ADM scheme proved to be numerically unstable, defining the 
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Table 3: A summary of groups reporting NS-NS merger calculation results. The asterisk for the 
KT collaboration's MHD column indicates that they have used an MHD-based code for other 
projects, but not yet for NS-NS merger simulations. Gravitational formalisms include full GR, 
assumed to be implemented using the BSSN decomposition except for the HAD collaborations's 
GHG approach, the CF approximation, or Newtonian gravity. Microphysical treatments include 
physically motivated EOS models or quark-matter EOS and neutrino leakage schemes. 
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GR 
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GR 
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CF 
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Quark, Phys. EOS, j/-leakage 


Bremen 


181^1 11^2911:^43] 


Newt 


Y 


Phys. EOS , i^-leakage 



auxiliary quantities f * = — 7*j and treating these expressions as independent variables stabilized the 
system and allowed for long-term evolutions. While slight variants exist, the 19 evolved variables 
are typically either the conformal factor ip or its logarithm (/), the conformal 3- metric jij , the trace 
K of the extrinsic curvature, the trace free extrinsic curvature Aij and the conformal connection 
functions P*. The evolution equations themselves are given in Appendix [XI 

The BSSN scheme was used in the binary merger calculations of the KT collaboration |2821 
12831 12801 1281] , the first completely successful NS-NS calculations ever performed in full GR. Ever 
since the UTB/RIT |61] and Goddard fS? groups showed simultaneously that a careful choice of 
gauge allows BHs to be evolved in BSSN schemes without the need to excise the singularity, these 
"puncture gauges" have gained widespread hold, and have been used to evolve NS-NS binaries 
(and in some cases, BH-NS binaries) by the KT coUaboration [317], UIUC [T7D], and Whisky p7] . 

The generalized harmonic formalism, developed over about two decades from initial theoretical 
suggestions up to its current numerical implementation |111| 1114) 1228| 1127| 1227] was used to 
perform the first calculations of merging BH-BH binaries by Pretorius [227], and has since been 
applied to NS-NS binaries by the HAD collaboration j7] and to BH-NS mergers by HAD [67] . 
SXS [57] IHH 1107] , and the Princeton group |2891 155] . It involves introducing a set of auxiliary 
quantities denoted /i'^ representing the action of the wave operator on the spacetime coordinates 
themselves 

i/^ = P^ = .g"^P^^ = __La,(y^.g'^'^) = .g^^V^V^x'^ ^ D (25) 

which are treated as independent gauge variables whose evolution equation must be specified. 
Current first-order formulations |1691 [7] evolve equations for the spacetime metric g^jy along with 
its spatial derivative, = digfj,„ and projected time derivative H^j^ = —n°'dag^v, subject to 

consistency constraints on the spatial derivatives. The first BH-NS merger calculations in the 
GH formalism used a first-order reduction |169] of the Einstein equations and specified the source 
functions to damp to zero exponentially in time, while the first binary NS merger work 7^ used a 
similar first-order reduction and chose harmonic coordinates with — 0. 

In both formalisms, most groups employ grid-based finite differencing to evaluate spatial deriva- 
tives. While finite differencing operators may be easily written down to arbitrary orders of accuracy, 
there is a trade-off between the computational efficiency achievable by using smaller, second-order 
stencils and the higher accuracy that can be attained using larger, higher-order ones. For the 
moment, many groups are now moving to at least fourth-order accurate differencing techniques, 
with a high likelihood that at least the field sector of NS-NS merger calculations will soon be 
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performed at comparable order to BH-BH calculations, at either sixth [135] or eighth-order |178j 
accuracy, if not higher. The main limitation to date involves the complexity of shock capturing 
using higher-order schemes, as we discuss in Section [5. 2.31 below. 

Imposing physically realistic and accurate boundary conditions remains a difficult task for 
numerical codes. Ideally, one wishes to impose boundary conditions at large distances that preserve 
the GR constraints and yield a well-posed initial-boundary value problem. On physical grounds, 
the boundary should only permit outgoing waves, preventing the reflection of spurious waves back 
into the numerical grid. Otherwise, reflections may be avoided by placing the outer boundaries 
so far away from the matter that they remain causally disconnected from the merging objects for 
the full duration of the simulation. Building upon previous work (see, e.g., |112| I150| [T2l 11491 
I252| 1238] . Winicour |324) derived boundary conditions that satisfy all of the desired conditions 
for the generalized harmonic formalism. No such conditions have been derived for the full BSSN 
formalism, though progress has been made (see, e.g., [i?lll28) ) so that we may now construct well- 
posed BCs in the weak-gravity linearized limit of BSSN [200] and for related first-order gravitational 
formalisms [49] . 

5.2.2 Grid decompositions 

While unigrid schemes are extremely convenient, they tend to be extremely inefficient, since one 
must resolve small-scale features in the central regions of a merger but also extend grids out to 
the point where the GW signal has taken on its roughly asymptotic form Thus, nearly every 
code incorporates some means of focusing resolution on the high-density regions via some form of 
mesh refinement. A simple approach, for instance, is to use "fisheye" coordinates that represent a 
continuous radial deformation of a grid |170) , a technique that had previously been used successfully, 
e.g., for BH-BH mergers pTIlS^ . 

While fixed mesh refinement offers the chance for greater computational efficiency and accu- 
racy, much current work focuses on adaptive mesh refinement, in which the level of refinement of 
the grid is allowed to evolve dynamically to react to the evolving binary configuration. Several 
codes, both public and private, now implement this functionality. The publicly available Carpet 
computational toolkit [2581 , which is distributed to the community as part of the open source 
Einstein Toolkit [SHI 1173] uses a "moving boxes" approach, and has been designed to be com- 
patible with the widely used and publicly available Cactus framework. It has been successfully 
implemented by the Whisky group [17] to perform NS~NS mergers, by UIUC for their BH-NS 
mergers [93], and a host of other groups for BH-BH mergers and additional problems. The KT 
code SACRA also implements an adaptive mesh refinement (AMR) scheme for NS-NS and BH-NS 
mergers [327] , as does the most recent version of the HAD collaboration's code [7] , which is based 
on the publicly available infrastructure of the same name |131l [6] , and the BAM code employed by 
the Jena group [3031 1411 1121] . The Princeton group also has an AMR code, which has been used 
to perform BH-NS mergers to date [2891 EHl M\ 

One drawback of employing rectangular grids is that memory costs scale like N^, where N is the 
number of grid cells across a side, and total computational effort like N'^ once Courant conditions 
are figured in. Since one does not necessarily need high angular resolution at large radii, there is 
great current interest in developing schemes that use some form of spheroidal grid, for which the 
memory scaling is merely oc N. A group at LSU has implemented a multi-patch method [332], in 
which space is broken up into a number of non-overlapping domains in such a way that the six 
outermost regions (projections of the faces of a cube onto spheres of constant radius), maintain 
constant angular resolution and thus produce linear dependence of the total number of grid points 
on the number of radial points. To date, it has been used primarily for vacuum spacetimes [220] 
and hydrodynamics on a fixed background. The SXS collaboration, begun at Caltech and Cornell 
and now including members at CITA and Washington State, has used a spectral evolution code 
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with multiple domains to evolve BH-NS binaries, which achieves the same scaling by expanding 
the metric fields in modes rather than in position space |HZ] • Their first published results on NS-NS 
binaries are currently in preparation (see [140] for a brief summary of work to date). 

While all of these grid techniques produce tremendous advantages in computational efficiency, 
each required careful study since deformations of a grid or the introduction of multiple domains 
can introduce inaccuracies and non-conservative effects. As an example, in AMR schemes, one 
must deal with the same reflection issues at reflnement boundaries that are present at the physical 
boundaries of the grid, as discussed above, though the interior nature of the boundaries allows for 
additional techniques, such as tapered grid boundaries [167j . to be used to minimize reflections 
there. The study of how to minimize spurious effects in these schemes continues, and will represent 
an important topic for years to come, especially as evolution schemes become more complicated 
by including more physical effects. 

5.2.3 Hydrodynamics, MHD, and high-resolution shock capturing 

Fluids cannot be treated in the same way as the spacetime metric because flnite differencing 
operators do not return meaningful results when placed across discontinuities induced by shocks. 
Instead, GR(M)HD calculations must include some form of shock-capturing that accounts for these 
jumps. These are typically implemented in conservative schemes, in which the fluid variables are 
transformed from the standard "primitive" set P, which includes the fluid density, pressure, and 
velocity (and magnetic fleld in MHD evolutions), into a new set U so that the evolution equations 
may be written in the form 

dt{U)^V -F + S, (26) 

where the flux functions F[P) and source terms S{P) can be expressed in terms of the primitive 
variables but not their derivatives. These schemes allow one to evolve the resulting MHD set of 
equations so that numerical fluxes are conserved to numerical precision across cell walls as the fluid 
evolves in time. One widely used scheme, often referred to as the Valencia formulation [184) . is 
described in Appendix IbI. 

There are important mathematical reasons for casting the GRHD/GRMHD system in conser- 
vative form, primarily since the mathematical techniques describing Godunov methods may be 
called into play |120) . In such methods, we assume that the evolution of the quantities U may be 
expressed in the form 

U{x = Xi,t^ tn+i) = U{x ^ Xt,t = t„) + ^ (^F{x = Xi_xli) - F{x = 2^4+1/2)) (27) 

where the points have spatial coordinates Xi = xa + iAx and the tiniesteps satisfy t„ = to + nAt. 
The fluxes must be determined by solving the Riemann problem at each cell face (thus the half- 
integer indices), either exactly or approximately. It can be shown that solutions constructed this 
way, when convergent, must converge to a solution of the original problem, even when shocks are 
present 162 . 

First one reconstructs the primitives from the conserved variables on both sides of an interface, 
using interpolation schemes designed to respect the presence of shocks. Simple schemes involving 
limiters yield second-order accuracy in general but first-order accuracy at shocks, while higher-order 
methods such as PPM (piecewise parabolic method) and essentially non-oscillatory (ENO) schemes 
such as CENG (central ENO) and WENO (weighted ENO) yield higher accuracy but at much 
higher computational cost. Once primitives are interpolated to the cell interfaces, flux terms are 
evaluated there and one solves the Riemann problem describing the evolution of two distinct fluid 
configurations placed on either side of a membrane (see [106] for a description). While complete 
solutions of the Riemann problem are painstaking to evolve, a number of approximation techniques 
exist and do not reduce the order of accuracy of the code. Finally, one must take the conservative 
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solution, now advanced forward in time, and recover the underlying primitive variables, a process 
that requires solving as many as eight simultaneous equations in the case of GRMHD or five for 
GRHD systems. A number of methods to do this have been widely studied [199) . and simplifying 
techniques are known for specific cases (for the case of polytropic EOSs in GRHD evolutions, one 
need only invert a single non-analytic expression and the remaining variables can then be derived 
analytically). 

The inclusion of magnetic fields in hydrodynamic calculations adds another layer of complex- 
ity beyond shock capturing. Magnetic fields must be evolved in such a way that they remain 
divergence- free, much in the same way that relativistic evolutions must satisfy the Hamiltonian and 
momentum constraints. Brute force attempts to subtract away any spurious divergence often lead 
to instabilities, so more intricate schemes have been developed. "Divergence cleaning" schemes typ- 
ically introduce a new field representing the magnetic field divergence and use parabolic/hyperbolic 
equations to damp the divergence away while moving it off the computational domain; the approach 
is relatively simple to implement but prone to small-scale numerical errors 24 . "Constrained 
transport" schemes stagger the grids on which different physical terms are calculated to enforce 
the constraints (see, e.g., [3051 for a particular implementation), and have been applied widely to 
many different physical configurations. Recently, a new scheme in which the vector potential is 
used rather than the magnetic field was introduced by Etienne, Liu, and Shapiro [921 1951 . and 
found to yield successful results for a variety of physical configurations including NSs and BHs. 



5.3 Microphysical numerical techniques 
5.3.1 Neutron star physics and equations of state 

One of the largest uncertainties in the input physics of NS-NS merger simulations is the true 
behavior of the nuclear matter EOS. To date, EM observations have yielded relatively weak con- 
straints on the NS mass-radius relation, with the most precise simultaneous measurement of both 
as of now resulting from observations of Type 1 X-ray bursts from accreting NSs in three different 
sources |216] . In each case, the NS mass was found to lie in the range 1.3 Mq < Mns ^ 2 Mq and 
the radius 8 km < i?NS ^ 12 km, implying a NS compactness 

C ^ = 0.1476 ( ^] ( ' « 0.16 - 0.37. (28) 

A more stringent constraint on the NS EOS is provided by observations of the Shapiro time 
delay in the binary millisecond pulsar PSR J1614-2230, which was found to have a mass Mns = 
1.97± 0.04 Mq [81], which would rule out extremely soft EOS models incapable of supporting such 
a massive NS against collapse. As we discuss in more detail below, GW observations are likely 
to eventually yield tighter constraints than our current EM-based ones, though BH-NS mergers, 
which can undergo stronger tidal disruptions than NS-NS mergers at frequencies closer to LIGO 
and other GW observatories' maximum frequency sensitivity band, may prove to be more useful 
for the task than NS-NS mergers. 

Given the large theoretical uncertainties in describing the proper physical NS EOS, many groups 
have chosen the simplest possible parameterization: a polytrope (see Eg. I12p . Under this choice, 
the enthalpy h takes the particularly simple form 

h=l + e + P/p^l^Te. (29) 

Initial data are generally assumed to follow the relation 

P = Kp^, (30) 
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where K is constant across the fluid. In the presence of shocks, the value of K for a particular fluid 
element will increase with time. We note that the Whisky group [T7] [TB] uses the term "polytropic" 
to refer to simulations in which Eg. 1301 is enforced throughout, which implies adiabatic evolution 
without shock heating, and use the term "ideal fluid" to describe an EOS that includes the effects 
of shock heating and enforces Eq. [29l 

Since the temperatures of NSs typically yield thermal energies per baryon substantially below 
the Fermi energy, one may treat nearly all NSs as effectively cold, except for the most recently 
born ones. During the merger process for NS-NS binaries, the matter will remain cold until the 
two NSs are tidally disrupted and a disk forms, at which point the thermal energy input and 
substantially reduced fluid densities require a temperature evolution model to properly model the 
underlying physics. In light of these results, some groups adopt a two-phase model for the NS EOS 
(see, e.g., j281j ). where a cold, zero-temperature EOS, evaluated as a function of the density only, 
encodes as much information about as we possess about the NS EOS, and the hot phase depends 
on both the density and internal energy, typically in a polytropic way 

F(p,£) = Pcold(p)+-Phot(p,£) [Phot(p,£) = (r- l)p(£-£cold)]. (31) 

There are a number of physically motivated EOS models that have been implemented for 
merger simulations, whose exact properties vary depending on the assumptions of the underlying 
model. These include models for which the pressure is tabulated as a function of the density 
only: EPS [218) . SLy [83], and APR [3j; as well as models including a temperature dependence: 
Shen [2631 1262) and Lattimer-Swesty |161j . A variety of models have been used to study the 
effects of quarks, kaons, and other condensates, which typically serve to soften the EOS, leading 
to reduced maximum masses and more compact NSs |219l I118| I226| I119| [23] E] . 

Given the variance among even the physically motivated EOS models, it has proven useful to 
parameterize known EOS models with a much more restricted set of parameters. In a series of 
works, a Milwaukee/Tokyo collaboration determined that essentially all current EOS models could 
be fit using four parameters, so that their imprint on GW signal properties could be easily analyzed 
|233l 12341 1182) . Their method assumes that the SLy EOS describes NS matter at low densities, 
and that the EOS at higher densities can be described by a piecewise polytropic fit with breaks at 
p = 10^^'^ and 10^^ g/cm'^. The four resulting parameters are Pi — P (^p — 10^^'^), the pressure 
at the first breakpoint density, which normalizes the overall density scale, as well as Fi, Fs, the 
adiabatic exponents in the three regions. Their results indicate that advanced LIGO should be 
able to determine the NS radius to approximately 1 km at an effective distance of 100 Mpc, which 
would place tight constraints on the value of Pi in particular. 

5.4 Electromagnetic and neutrino signature modeling 

Motivated by the evidence that SGRBs frequently appear in galaxies with very low star forma- 
tion rates [40] 1109) , astronomers have suggested that their progenitors are likely to be mergers of 
either NS-NS and/or BH-NS binaries. While soft-gamma repeaters (SGRs) have been confirmed 
as an SGRB source from observations of the system SGR 1806-20, they make up no more than 
approximately 15% of the total observed SGRB fraction according to the leading population es- 
timates |163l 1195) . There has been much interest in predicting the EM signatures of NS-NS and 
BH-NS mergers, along with the associated neutrino emission. The simplest models estimate a 
local radiation cooling rate for the matter but do not attempt to follow the paths of the photons 
and/or neutrinos after they are emitted, instead calculating the time-dependent luminosity assum- 
ing free streaming. Such models have been used in non-GR simulations of binary mergers going 
back more than a decade |249l 1242) , and recently such schemes have been used to perform full GR 
NS-NS mergers [260] , including a self-consistent evolution of the electron fraction of the material 
Ye, rather than a passive advection approach. 
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More complicated flux-limited diffusion schemes, in which the neutrino fluxes for given species 
and energies are given by explicit formulae that limit to the correct values for zero optical depth 
(free-streaming) and very large optical depth (diffusion), have been used as a post-processing tool 
to investigate the merger remnants in Newtonian NS-NS mergers [82], but have yet to be applied 
to full GR simulations. Finally, radiation transport schemes to evolve EM and neutrino fluxes 
passing through fluid configurations have been implemented in numerical GR codes |80[ 1103] , but 
have yet to be used in binary merger simulations. 



5.5 GW signal modeling 

Measuring the GW signal from a dynamical merger calculation is a rather difficult task. One must 
determine, using a method unaffected by gauge effects, the perturbations at asymptotically large 
distances from a source by extrapolating various quantities measured at large but finite distances 
from the merger itself. 

In the early days of numerical merger simulations, most groups typically assumed Newtonian 
and/or quasi- Newtonian gravitation, for which there is no well-defined dynamical spacetime met- 
ric. GW signals were typically calculated using the quadrupole formalism, which technically only 
applies for slow-moving, non-relativistic sources (see [189] for a thorough review of the theory). 
Temporarily reintroducing physical constants, the strains of the two polarizations for signals emit- 
ted in the z-direction are 

G 

h - 

where r is the distance from the source to the observer and i^- it the traceless quadrupole moment 
of the system. The GW luminosity and angular momentum loss rate of the system are given, 
respectively, by 

'^^ ) GW 

dLk \ _ 2G ■■■■■ 

While only approximate, the quadrupole formulae do yield equations that are extremely straight- 
forward to implement in both grid and particle-based codes using standard integration techniques. 

Quadrupole methods were adopted for later PN and CF simulations, again because the metric 
was assumed either to be static or artificially constrained in such a way that made self-consistent 
determination of the GW signal impossible. One important development from this period was the 
introduction of a simple method to calculate the GW energy spectrum dE/df from the GW time- 
series through Fourier transforming into the frequency domain [325]. GW signals analyzed in the 
frequency domain allowed for direct comparison with the LIGO noise curve, making it much easier 
to determine approximate distances at which various GW sources would be detectable and the 
potential signal-to-noise ratio that would result from a template search. To constrain the nuclear 
matter EOS, one can examine where a GW merger spectrum deviates in a measurable way from 
the quadrupole point-mass form, 

\~^] = 5 (7rG(TOi -I- 1712) Jgw) ; /gw = 2/orb = , (32) 

\df Jgw ^ ^ 

because of finite-size effects, and then link the deviation to the properties of the NS 'SB], as we 
show in Figure [9] 
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Figure 9: Approximate energy spectrum dEcw/df derived from QE sequences of equal- mass NS- 
NS binaries with isolated ADM masses Mns = 1.35 and a F = 2 EOS, but varying com- 
pactnesses (denoted M/R here), originally described in 2981. The diagonal lines show the energy 
spectrum corresponding to a point-mass binary, as well as values with 90%, 75%, and 50% of the 
power at a given frequency (see Figure 2 of [98i)- Asterisks indicate the onset of mass-shedding, 
beyond which QE results are no longer valid. 



Full GR dynamical calculations, in which the metric is evolved according to the Einstein equa- 
tions, generally use one of two approaches to calculate the GW signal from the merger, if not both. 
The first method, developed first by by Regge and Wheeler |235| and Zerilli |330J and written 
down in a gauge-invariant way by Moncrief |191J involves analyzing perturbations of the metric 
away from a Schwarzschild background. The second uses the Newman-Penrose formalism [198] to 
calculate the Weyl scalar tlj4, a contraction of the Weyl curvature tensor, to represent the outgoing 
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wave content on a specially constructed null tetrad that may be calculated approximately |60j . 
The two methods are complementary since they incorporate different metric information and re- 
quire different numerical integrations to produce a GW time series. Regardless of the method 
used to calculate the GW signal, results are often presented by calculating the dominant s = —2 
spin- weighted spherical harmonic mode. For circular binaries, the I = 2, m = 2 mode generally 
carries the most energy, followed by other harmonics; in cases where the components of the binary 
have nearly equal masses and the orbit is circular, the falloff is typically quite rapid, while extreme 
mass ratios can pump a significant amount of the total energy into other harmonics. For elliptical 
orbits, other modes can dominate the signal, e.g., a 3 : 1 ratio in power for the 1 = 2, m = mode 
to the I = 2, m = 2 mode observed for high-ellipticity close orbits in [121) . A thorough summary 
of both methods and their implementation may be found in [253J . 
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6 Dynamical Calculations 



NS~NS merger simulations address a broad set of questions, which can be roughly summarized as 
follows (note the same questions apply to BH-NS mergers as well): 

1. What is the final fate of the system, assuming a given set of initial parameters? Do we 
get a prompt collapse to a BH or the formation of a HMNS supported against collapse by 
differential rotation? Other outcomes are disfavored, at least for pre-merger NSs with masses 
Mns > 1.4 Mq since the supramassive limit is at most 20% larger than that of a non-rotating 
NS [7T][72], and even for the stiffest EOS these values are typically less than 2.8 A/q. 

2. What is the GW signal from the merger, and how does it inform us about the initial pre- 
merger parameters of the system? 

3. What fraction of the system mass is left in a disk around the central BH or HMNS? While 
deriving exact EM emission profiles from a hydrodynamical configuration remains a challenge 
for the future, minimum conditions that would allow for the energy release observed in SGRBs 
have been established based on scaling arguments. 

4. What is the neutrino and EM emission from the system, in both the time and energy domains? 
Obviously, the answer to this question and those that follow depend critically on the answers 
above. 

5. What role do B-fields play in the GW, EM, and neutrino emission, and how does that tie in 
with other models suspected of having the same disk/jet geometry and gamma-ray emission 
like active galactic nuclei, pre-main sequence stars, etc.? 

6. Do mergers produce a cosmologically significant quantity of r-process elements, or do those 
likely get produced by other astrophysical events instead? 

The influence of the gravitational formalism used in a numerical simulation on the answer one 
finds for the questions above differs item by item. Determining the final fate of a merging system 
is highly dependent on the gravitational formalism; NS-NS merger remnants only undergo collapse 
in quasi-relativistic and fully GR schemes. Moreover, orbital dynamics at separations comparable 
to the ISCO and even somewhat larger depend strongly on the gravitational scheme. In particular, 
mass loss rates into a disk are often suppressed by orders of magnitude in GR calculations when 
compared to CF simulations, and even more so in comparison to PN and Newtonian calculations. 
EM emission profiles from a disk are difficult to calculate accurately without the use of full GR 
for this reason. On the other hand, while GR is required to calculate the exact GW signal from a 
merger, even early Newtonian simulations predicted many of the qualitative GW emission features 
correctly, and PN and CF schemes yielded results with some degree of quantitative accuracy about 
the full wavetrain. B-fields have only begun to be explored, but it already seems clear that they will 
affect the hydrodynamical evolution primarily after the merger in cases where differential rotation 
in a HMNS or disk winds up magnetic field strengths up to energy equipartition levels, vastly 
stronger than those found in pre-merger NSs. For such configurations, non-relativistic calculations 
can often reproduce the basic physical scenario but full GR is required to properly understand 
the underlying dynamics. Finally, the production of r-process elements, which depends sensitively 
on the thermodynamic evolution of the merger, seems to generally disfavor binary mergers as a 
significant source of the observed stellar abundances since the temperature and thus the electron 
fraction of the fiuid remains too small |248| , regardless of the nature of the gravitational treatment 
used in the calculations. This picture may need to be revised if significant mass loss occurs from 
the hot accretion disk that forms around the central post- merger object, but numerical calculations 
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do not currently predict sufficient mass loss to match observations |292j . We will address each of 
these topics in greater detail in the sections below. 

Since the first NS-NS merger calculations, there have been two main directions for improve- 
ments: more accurate relativistic gravitation, resulting in the current codes that operate using a 
self-consistent fully GR approach, and the addition of microphysical effects, which now include 
treatments of magnetic fields and neutrino/EM radiation. Noting that several of the following 
developments overlapped in time, e.g., the first full GR simulations by Shibata and Uryii |282j 
are coincident with the first PN SPH calculations, and predate the first CF SPH calculations, we 
consider in turn the original Newtonian calculations, those performed using approximate relativis- 
tic schemes, the calculations performed using full GR, and finally those that have included more 
advanced microphysical treatments. 

6.1 Quasi-equilibrium and semi-analytic methods vs fully dynamical re- 
sults 

Before reviewing fully dynamical calculations of NS-NS mergers, it is worthwhile to ask how much 
information can already be deduced from QE calculations, which may be performed at much 
smaller computational cost, as well as from semi-analytic PN treatments and related approximate 
techniques. Clearly, the details of the merger and ringdown phases fall outside the QE regime, 
so only dynamical calculations can yield reliable information about the stability of remnants, 
properties of ejecta, or other processes that arise during the merger itself or in its aftermath. 
Thus, the primary point of comparison is the GW signal just prior to merger, which is also easier 
to detect (for first and second generation interferometers). 

The strength of QE calculations lies in their ability to model self-consistently finite-size effects 
not captured in PN treatments (which always assume two orbiting point masses). The increased 
tidal interaction between the objects typically results in a more rapid phase advance of the binary 
orbit, which is important for constructing template waveforms that cover the entire NS-NS inspiral, 
merger, and ringdown. While QE sequences potentially offer a wealth of information about well- 
separated binaries and can help fix the phase evolution of the inspiraling binary, they do have two 
weaknesses arising as the binary approaches the stability limit. First, most QE methods, including 
the GTS formalism described in Sec. 14.2.11 are time-symmetric, and assume that the NS possess a 
symmetry plane perpendicular to the direction of motion (i.e., a front-back symmetry whose axis 
is perpendicular to the orbital angular momentum and the binary separation vector). In reality, 
tidal lags develop prior to final plunge, with the innermost edge of each NS rotating forward and 
the outer edge backwards. This effect has been captured in analytic and semi-analytic approaches 
(see, e.g., 1158' for an early example), and is clearly seen in dynamical calculations (see Fig. [3]), but 
is not captured in GTS-based schemes (tidal lags also develop in BH~NS merger calculations when 
the BH has a non-zero spin, since this breaks the front-back symmetry; see [297J for an example). 

A second weakness of QE methods is the treatment of the ISCO, particularly its importance 
as a characteristic point along an evolutionary sequence that, in theory, could encode information 
about the NS EOS. Simple estimates of the infall trajectory derived solely from QE sequences 
predict a very sudden and rapid infall near the ISGO, i.e., the point where the binding energy 
reaches a minimum along the sequence (see, e.g., the argument in [98]). However, this is clearly 
an oversimplification. In reality, binaries transition more gradually to the merger phase, and the 
inward plunge may occur significantly before reaching the formal ISGO; this in turns leads to more 
rapidly growing deviations from the QE approximation. Looking at the GW energy spectrum, 
one typically sees minor deviations from the point-mass predictions at frequencies below those 
characterizing the ISGO, but substantially more power at frequencies above it. Equivalently, the 
cutoff frequency for GW emission /cut, where the spectrum starts deviating strongly from the 
point- mass prediction, is usually higher than the QE frequency near the ISGO, /iscOj while simple 
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QE estimates assume these two frequencies to coincide. 

To date, most attempts to generate waveforms extended back to arbitrarily early starting points 
involve numerically matching PN signals, typically generated using the Taylor T4 approach [53] . 
onto the early stages of numerically generated waveforms, with some form of maximum overlap 
method used to provide the most physical transition from one to the other. These approaches may 
be improved by adding tidal effects to the evolution, typically parameterized by the tidal Love 
numbers that describe how tidal gravity fields induce quadrupole deformations |105j . Tidal effects 
can be placed into a relativistic framework [551 El] ) which may be included within the effective 
one-body (EOB) formalism to produce high-accuracy waveforms [75|. In the EOB approach [58] , 
resummation methods are used to include higher-order PN effects, though some otherwise unfixed 
parameters need to be set by comparing to numerical simulations. 

Work is in its early stages to compare directly the GW spectra inferred from QE sequences 
of NS-NS binaries with those generated in numerical relativity simulations, but this comparison 
has been discussed at some length with regard to BH-NS mergers. Noting that NS-NS mergers 
generally correspond more closely to the BH-NS cases in which an ISCO is reached prior to the 
onset of tidal disruption, the KT collaboration |278| I271j concluded that the cutoff frequency 
marking significant deviations from PN point-mass behavior is roughly 30% higher than that 
marking emission near the classical ISCO for BH-NS systems (/cut — 1-3/isco)- 

A more detailed study has now been performed comparing EOB methods to numerical evo- 
lutions. By comparing to long-term simulations of NS-NS mergers, Baiotti et al. [15 find that 
EOB models may be tuned, via careful choices of their unfixed parameters, to reproduce the GW 
phases and amplitudes seen in NR evolutions up until the onset of merger. They further suggest 
that the EOB approach seems to cover a wider range of phase space than the Taylor T4 approach, 
presumable because of a more consistent representation of tidal effects, and offers the best route 
forward for construction of more accurate NS-NS inspiral templates. 

6.2 Early dynamical calculations 

The earliest NS-NS merger calculations were performed in Newtonian gravity, sometimes with the 
addition of lowest-order 2.5PN radiation reaction forces, and typically assumed that the NS EOS 
was polytropic. Both Eulerian grid codes [210, 131, 2IlS [IMl EISl ESfl IWl |293] and Lagrangian 
SPH codes I23I1 12311 Hi MB were employed, and GW signals were derived under the 
assumptions of the quadrupole formalism. Configurations in Newtonian gravity cannot collapse, 
so a stable (possibly hypermassive) remnant was always formed. For polytropic EOS models with 
adiabatic indices larger than the classical minimum for production of a Jacobi ellipsoid, F > 2.6, 
remnants were typically triaxial and maintained a significant-amplitude GW signal until the end of 
the simulation. For simulations using smaller values of F, remnants rapidly relaxed to spheroidal 
configurations, quickly damping away the resulting GW signal. Mass loss from the central remnant 
was often quite significant, with thick accretion disks or completely unbound material comprising 
up to 10-20% of the total system mass 

Mass loss was suppressed in numerical simulations by constructing irrotational, rather than 
synchronized, initial data. Irrotational flow is widely thought to be the more physically realistic 
case, since viscous forces are much too weak to synchronize a NS prior to merger [451 1145] . When 
irrotational NSs (which are counter-rotating in the corotating frame of the binary) first make con- 
tact, a vortex sheet forms. Since the low-density fluid layers at the contact surface are surrounded 
at first contact by the denser fluid layers located originally within each NS, the conflguration is well 
understood to be Kelvin-Helmholtz unstable, resulting in rapid mixing through vortex production. 
Meanwhile, mass loss through the outer Lagrange points is hampered by the reduced rotational 
velocity along the outer halves of each NS. 

The GW emission from these mergers is composed of a "chirp," increasing in frequency and 
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amplitude as the NSs spiral inward, followed by a ringdown signal once the stars collide and merge. 
In |325ll326j . a procedure to calculate the energy spectrum in the frequency band was laid out, with 
the resulting signal following the quadrupole, point-mass power-law form up to GW frequencies 
characterizing the beginning of the plunge. Above the plunge frequency, a sharp drop in the GW 
energy was seen, followed in some cases by spikes at kHz frequencies representing coherent emission 
during the ringdown phase. 

6.3 Approximate relativistic schemes 

The first steps toward approximating the effects of GR included the use of IPN dynamics or the 
CF approximation. Using a formalism derived by Blanchet, Damour, and Schaefer [35], the IPN 
equations of motion require the solution of eight Poisson-like equations in the form 

VV = /(a?), (33) 

where the source terms f{x) are compactly supported, and thus the fields "0 may be determined 
using the same techniques already in place to find the Newtonian potential. Adding in the 
lowest-order dissipative radiation reaction effects requires solution of a ninth Poisson equation 
for a reaction potential. The IPN formalism was implemented in both grid-based [212] and SPH 
codes [lOl [99] IIOII 1100] . Unfortunately, physically realistic NSs are difficult to model using a PN 
expansion, since the characteristic NS compactness C > 0.15, leads to first order "corrections" that 
often rival Newtonian terms in magnitude. To deal with this problem Ayal et al. TO^ considered 
large (R 30 km), low-mass (< 1 Mq) NSs, allowing them to study relativistic effects but making 
results more difficult to interpret for physically realistic mergers. In [99l IIOII 1100) . a dual speed of 
light approach was used, in which all IPN effects were scaled down by a constant factor to yield 
smaller quantities while Newtonian and radiation reaction terms were included at full-strength. 
Both SPH groups found that the GW signal in PN mergers is strongly modulated, whereas New- 
tonian merger calculations typically yielded smooth, either monotonically decreasing or nearly 
constant-amplitude ringdown signals. Even reduced IPN effects were shown to suppress mass loss 
by a factor of 2 - 5 for initially synchronized cases, and disk formation was seen to be virtually 
non-existent for initially irrotational, equal-mass NSs with a stiff (P = 3 polytropic) EOS [101) . 

Moving beyond the linearized regime, several groups explored the CF approximation, which 
incorporates many of the nonlinear effects of GR into an elliptic, rather than hyperbolic, evolution 
scheme. While non-linear elliptic solvers are expensive computationally, they typically yield stable 
evolution schemes since field solutions are always calculated instantaneously from the given matter 
configuration. Summarized quickly, the CF approach involves solving the CTS field equations, 
Eqs. 1181 1191 and 1201 ^-t every timestep, and evolving the matter configuration forward in time. 
The metric fields act like potentials, with various gradients appearing in the Euler and energy 
equations. While the CTS formalism remains the most widely used method to construct NS-NS 
(and BH-NS) initial data, it does not provide a completely consistent dynamical solution to the 
GR field equations. In particular, while it reproduces spherically symmetric configurations like 
the Schwarzschild solution exactly, it cannot describe more complicated configurations, including 
Kerr BHs. Moreover, because the CF approximation is time-symmetric, it also does not allow one 
to consistently predict the GW signal from a merging configuration. As a result, most dynamical 
calculations are performed by adding the lowest-order dissipative radiation reaction terms, either 
in the quadrupole limit or via the radiation reaction potential introduced in [48] . 

The CTS equations themselves were originally written down in essentially complete form by 
Isenberg in the 1970s, but his paper was rejected and only published after a delay of nearly 30 
years |136] . In the intervening years, Wilson, Mathews, and Marronetti [3221 13231 11861 [T55] inde- 
pendently re-derived the entire formalism and used it to perform the first non-linear calculations 
of NS-NS mergers (as a result, the formalism is often referred to as the "Wilson-Mathews" or 
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"Isenberg- Wilson-Mathews" formalism) . The key result in j322l I323| 11861 1185J was the existence 
of a "collapse instability," in which the deeper gravitational wells experienced by the NSs as they 
approached each other prior to merger could force one or both to collapse to BHs prior to the 
orbit itself becoming unstable. Unfortunately, their results were affected by an error, pointed out 
in [104) . which meant that much of the observed compression was spurious. While their later cal- 
culations still found some increase in the central density as the NSs approached each other [187] . 
these results have been contradicted by other QE sequence calculations (see, e.g., [308J ). Further- 
more, using a "CF-like" formalism in which the non-linear source terms for the field equations 
are ignored, dynamical calculations demonstrated the maximum allowed mass for a NS actually 
increases in response to the growing tidal stress [268) . 

The CP approach was adapted into a Lagrangian scheme for SPH calculations by the same 
groups that had investigated PN NS-NS mergers, with Oechslin, Rosswog, and Thielemann [207] 
using a multigrid scheme and Paber, Grandclement, and Rasio [97] a spectral solver based on the 
Lorene libraries [176j . The effects of non- linear gravity were immediately evident in both sets of 
calculations. In [207] . NS-NS binaries consisting of initially synchronized NSs merged without ap- 
preciable mass loss, with no more that ~ 10"'* of the total system mass ejected, strikingly different 
from previous Newtonian and PN simulations. When evolving initially irrotational systems, [97] 
found no appreciable developments of "spiral arms" whatsoever, indicating a complete lack of mass 
loss through the outer Lagrange points. Both groups also found strong emission from remnants 
for a stiff EOS, as the triaxial merger remnant produced an extended period of strong ringdown 
emission. Neither set of calculations indicated that the remnant should collapse promptly to form 
a BH, but given the high spin of the remnant it was noted that conformal flatness would have 
already broken down for those systems. 

6.4 Full GR calculations 

A summary of full GR calculations of NS-NS mergers is presented in Table 16.41 The KT col- 
laboration was responsible for the only full GR calculations of NS-NS mergers that predate the 
breakthrough calculations of numerically stable binary BH evolutions [2271 [22] IHT] . which have 
since transformed the field of GR hydrodynamics and MHD in addition to vacuum relativistic 
evolutions (Miller et al. [188] performed NS-NS inspiral calculations in full GR, but were not able 
to follow binaries through to merger). The first calculations of NS-NS mergers using a completely 
self-consistent treatment of GR were performed by Shibata and collaborators in the KT collab- 
oration using a grid based code and the BSSN formalism ,282j . GTS initial data consisting of 
equal-mass NSs described by a P = 2 polytropic EOS were constructed via SGP techniques [308] . 
for both synchronized and irrotational configurations. The hyperbolic system was evolved on a 
grid, with an approximate maximal slicing condition that results in a parabolic equation for the 
lapse [266] and an approximate minimal distortion condition for the shift vector requiring the solu- 
tion of an elliptic equation at every time step [267] . The shift vector gauge condition was found to 
fail when BHs were produced in the merger remnant, a well-known problem that had long bedev- 
iled simulations involving binary and even single BH evolutions, so modifications were introduced 
to extend the stability of the algorithm as far as possible. Among the key results from this early 
work was a clear differentiation between mergers of moderately low-compactness NSs (C > 0.11), 
where the remnant collapsed promptly to a BH, and very low-compactness models, which yielded 
hypermassive remnants stabilized against gravitational collapse by differential rotation. Virtually 
all the NS matter was contained within the remnant for initially irrotational models, which served 
as evidence against equal-mass NSs mergers being a leading source of r-process elements in the 
universe through ejection. The lack of significant mass loss in equal-mass mergers, together with 
insignificant shock-heating of the material, also argued against the likelihood of such mergers as 
progenitors for SGRBs if the gamma-ray emission was assumed to be coincident with the GW 
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Table 4: A summary of Full GR NS~NS merger calculations. EOS models include polytropes, piece- 
wise polytropes (PP), as well as physically motivated models including cold SLy \83^, FPS [218) . and 
APR [3] models to which one adds an ideal-gas hot component to reflect shock heating, as well as 
the Shen [263| 1262] finite temperature model and EOS that include Hyperonic contributions [259] . 
"Co/Ir" indicates that both corotating and irrotational models were considered; "BHB" indicates 
that BH binary mergers were also presented, including both BH-BH and BH-NS types, 'V-leak" 
indicates a neutrino leakage scheme was included in the calculation, "GH" indicates calculations 
were performed using the GHG formalism rather than BSSN, "non-QE" indicates superposition 
initial data were used, including cases where eccentric configurations were studied ("Eccen."); 
"MHD" indicates MHD was used to evolve the system. 



Group 


Ref. 


NS EOS 


Mass ratio 


C 


notes 


KT 


[282] 


F = 2 


1 


0.09-0.15 


Co/Ir 








F = 2,2.25 


0.89-1 


0.1-0.17 








|280| 


F = 2 


0.85-1 


0.1-0.12 








[281] 


SLy,FPS+Hot 


0.92-1 


0.1-0.13 








|277| 


SLy,APR-t-Hot 


0.64-1 


0.11-0.13 








|327| 


F = 2 


0.85-1 


0.14-0.16 


BHB 






|143| 


APR-t-Hot 


0.8-1 


0.14-0.18 








llMl 


APR,SLy,FPS+Hot 


0.8-1.0 


0.16-0.2 








|260| 


Shen 


1 


0.14-0.16 


i^-leak 






1133J 


PP+hot 


1 


0.12-0.17 










Shen, Hyp 


1.0 


0.14-0.16 


z/-leak 




HAD 


m 


F = 2 


1.0 


0.08 


GH, non- 


QE 




m 


F = 2 


1.0 


0.08 


GH, non- 


QE, MHD 


Whisky 


m 


F = 2 


1.0 


0.14-0.18 








m 


F = 2 


1.0 


0.20 








|115| 


F = 2 


1.0 


0.14-0.18 


MHD 






m 


F = 2 


1.0 


0.14-0.18 


MHD 








F = 2 


0.70-1.0 


0.09-0.17 








[a US] 


F = 2 


1.0 


0.12-0.14 








[237] 


F = 2 


1.0 


0.18 


MHD 




UIUC 


m 


F = 2 


0.85-1 


0.14-0.18 


MHD 


Jena 


[30311411 


F = 2 


1.0 


0.14 








[121] 


F = 2 


1.0 


1.4 


Eccen. 





39 



burst; instead a delayed burst following the collapse of a HMNS to a BH appeared more likely. 

Later works, in particular a paper by Shibata, Taniguchi, and Uryii |280) . introduced several 
new techniques to perform dynamical calculations that most codes at present still include in nearly 
the same or lightly modified form. These included the use of a high-resolution shock-capturing 
scheme for the hydrodynamics, as well as a Gamma-driver shift condition closely resembling the 
moving puncture gauge conditions that later proved instrumental in allowing for long-term BH 
evolution calculations. In the series of papers that followed their original calculations, the KT 
group established a number of results about NS-NS mergers that form the basis for much of our 
thinking about their hydrodynamic evolution: 

• By varying the EOS model for the NS as well as the mass ratio, it was possible to constrain 
the binary parameters separating cases that form a HMNS rather than producing prompt 
collapse to a BH, and it was quickly determined that the total system mass as a proportion 
of the maximum allowed mass for an isolated NS is the key parameter, with only weak 
dependence on the binary mass ratio. 

For polytropic EOS models, the critical compactness values leading to prompt collapse for 
equal-mass binary mergers were found to be C = 0.14 for F = 2 and C = 0.16 for F = 2.25 
[283] . As a rough rule, collapse occurred for polytropic EOS when the total system rest-mass 
was at least l.TMmax, where Mmax is the maximum mass of an isolated non-rotating NS for 
the given EOS. For physically motivated EOS models |281j . the critical mass was significantly 
smaller; indeed, the critical NS mass was found to be ~ 1.35 Mmax for the SLy EOS '83] (i.e., 
collapse for Mft > 2.7 Mq with M^ax = 2.04 Mq) and ~ 1.39 M^ax for the EPS EOS [2T8] 
(collapse for Mtot > 2.5 with Mmax — 1.8 M©). This was not a complete surprise, since 
for the physically motivated EOS the NS radius is nearly independent of the mass across 
much of the parameter space, limiting the ability of the HMNS to expand in response to the 
extra mass absorbed during the merger. 

• The mass ratio was found to play a critical role in the evolution of the remnant / disk con- 
figuration, since unequal-mass cases are better characterized as disruptions of the smaller 
secondary followed by its accretion onto the primary, rather than a true merger between 
the two NSs. Disk masses from full GR calculations |280| are generically smaller than those 
predicted from non-GR calculations. For polytropic EOS, disks contain approximately 4% of 
the total system mass for mass ratios q ~ 0.8, varying roughly cx (1 — g) for a fixed total mass, 
with the disk mass decreasing for heavier binaries (and thus larger compactnesses) given the 
stronger gravity of the central remnant. Using the stiffer APR EOS 0], the dependence on 
the mass ratio was seen to be much steeper for a physical EOS than for polytropes, scaling 
like (1 - qY, where p ~ 3 - 4 p77] . 

• With respect to GW emission, it was determined in |283) that in low-mass cases in which a 
HMNS was formed, stiffer polytropic EOS models were able to support the development of a 
bar-mode instability, leading to transient spiral arm formation from the remnant and an ex- 
tended period of strong GW emission, in the characteristic modulated form that results from 
differentially rotating ellipsoids (see, e.g., [159| I269| I255| [16]). Unequal-mass cases typically 
yielded one high frequency peak at roughly /gw ~ 2kHz corresponding to non-axisymmetric 
oscillations, and equal-mass cases yielded multiple peaks including those associated with 
quasi-radial oscillations as well [280]. For physical EOS models |281| . mass loss into a disk 
is reduced relative to the polytropic case given the higher compactness of the central region, 
and GW oscillation peaks, while very strong, occur at correspondingly higher frequencies. 
By contrast, prompt formation of a BH led to a ringdown signal with rapidly decreasing 
amplitude becoming negligible within a few dynamical times. 
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The GW signals were evaluated under the gauge-dependent assumption of transverse trace- 
lessness, and energy and angular momentum loss rates into each spherical harmonic mode 
were computed using the gauge-invariant Zerilli-Moncrief formalism |235l 13301 1191) in much 
the same way that is used by some groups in numerical relativity today (many BH-BH and 
hydrodynamics simulations report GW signals derived from the alternate -04 Weyl scalar 
formulation |198| l60] , or use both methods) . 

• In |281| . it was concluded that mergers of NSs with comparable masses made poor SGRB 
progenitor candidates, assuming prompt emission (because of the lack of energy available for 
neutrino annihilation), but that the energy budget in the HMNS case is orders of magnitude 
larger. Remarkably, this discussion from 2005 predates the first identifications of SGRBs with 
older populations, which greatly improved our theoretical understanding of compact object 
mergers as their likely progenitors. In the first work that followed the initial localizations of 
SGRBs, mergers of binaries with relatively small mass ratios, q 0.7, were seen to form suf- 
ficiently hot and massive disk to power a SGRB, albeit a relatively brief, low-luminosity one. 
It was suggested that the more likely SGRB progenitor is indeed a HMNS, since dissipative 
effects within the remnant can boost temperatures up to ~ 10^^ K. 

Further approximate relativistic investigations of NS-NS mergers, along with BH-NS merg- 
ers, as potential SGRB sources quickly swept through the community after the initial local- 
izations of SGRBs, with several groups using a wide variety of methods all concluding that 
mergers were plausible progenitors, but finding it extremely difficult to constrain the scenario 
in quantitative ways given the extremely complicated microphysics ultimately responsible for 
powering the burst (see, e.g., [203| I202j who investigated potential disk energies; [,239^, who 
modeled the fall-back accretion phase onto a BH; and |261', who considered the thermody- 
namic and nuclear evolution of disks around newly- formed BHs produced by mergers). We 
will return to this topic below in light of recent GRMHD Simulations. 

In the past few years, five groups have reported results from NS-NS mergers in full GR; KT, 
HAD, Whisky, UIUC, and Jena. Much of the work of the HAD and Whisky groups, developers 
respectively of the code of those names, began at Louisiana State University (HAD) and the Albert- 
Einstein Institut in Potsdam (Whisky), though both efforts now include several other collaborating 
institutions. Two other groups, the SXS collaboration that originated at Caltech and Cornell, 
and the Princeton group, have reported BH-NS merger results and are actively studying NS-NS 
mergers as well, but have yet to publish their initial papers about the latter. All of the current 
groups use AMR-based Eulerian grid codes, with four evolving Einstein's equations using the BSSN 
formalism and the HAD collaboration making use of the GHG method instead. HAD, Whisky, and 
UIUC have all reported results about magnetized NS-NS mergers (the KT collaboration has used 
a GRMHD code to study the evolution of magnetized HMNS, but not complete NS-NS mergers). 
The KT collaboration has considered a wide range of EOS models, including finite-temperature 
physical models such as the Shen EOS, and have also implemented a neutrino leakage scheme, 
while all other results reported to date have assumed a F = 2 polytropic EOS model. 

Given the similarities of the various codes used to study NS-NS mergers, it is worthwhile to 
ask whether they do produce consistent results. A comparison paper between the Whisky code 
and the KT collaboration's SACRA codes [50] found that both codes performed well for conservative 
global quantities, with global extrema such as the maximum rest-mass density in agreement to 
within 1% and waveform amplitudes and frequencies differing by no more than 10% throughout a 
full simulation, and typically much less. 

Several of the the groups listed above have also been leaders in the field of BH-NS simulations: 
the KT, HAD, and UIUC groups have all presented BH-NS merger results, as have the SXS 
collaboration [571 IMl [TUT] . and Princeton group PSU','55] (see fI7^ for a thorough review). 

We discuss the current understanding of NS-NS mergers in light of all these calculations below. 
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6.4.1 HMNS and BH remnant properties 



Using their newly developed SACRA code |327| , the KT group [143. , found that when a hybrid EOS 
is used to model the NS, in which the cold part is described by the APR EOS and the thermal 
component as a F = 2 ideal gas, the critical total binary mass for prompt collapse to a BH is 
-^tot = 2.8 — 2.9 Mq, independent of the initial binary mass ratio, a result consistent with previous 
explorations of other polytropic and physically motivated NS EOS models (see above). In all cases, 
the BH was formed with a spin parameter a w 0.78 depending very weakly on the total system 
mass and mass ratio. 

They further classified the critical masses for a number of other physical EOS in |133j . finding 
that binaries with total masses Aftot ^ 2.7 M0 should yield long-lived HMNSs (> 10 ms) and 
substantial disk masses with Mdisk > 0.04 M0 assuming that the current limit on the heaviest 
observed NS, M — 1.97 Mq [5T] is correct. In Figure [TU1 we show the final fate of the merger 
remnant as a function of the total pre-merger mass of the binary. "Type I" indicates a prompt 
collapse of the merger remnant to a BH, "Type H" a short-lived HMNS, which lasts for less than 
5 ms after the merger imtil its collapse, and "Type HI" a long-lived HMNS which survives for at 
least 5 ms. See [133] for an explanation of the EOS used in each simulation. 

3.0 '* 



2.9 



2.8 



2.7 



APR4 SLy ALF2 H3 H4 PS 
(11.428) (11.736) (13.188) (13.840) (13.759) (15.472) 

Figure 10: Type of final remnant corresponding to different EOS models (see Figure 3 of .133, ). 
The vertical axis shows the total mass of two NSs. The horizontal axis shows the EOSs together 
with the corresponding NS radii for A/ns — IAMq. 

While all of the above results incorporated shock heating, the addition of both finite-temperature 
effects in the EOS and neutrino emission modifies the numerically determined critical masses sep- 
arating HMNS formation from prompt collapse. Adding in a neutrino leakage scheme for a NS-NS 
merger performed using the relatively stiff finite-temperature Shen EOS, the KT collaboration re- 
ports in [260' that HMNSs will form generically for binary masses < 3.2 M©, not because they are 
centrifugally supported but rather because they are pressure-supported, with a remnant tempera- 
ture in the range 30 - 70 MeV. Since they are not supported by differential rotation, these HMNSs 
were predicted to be stable until neutrino cooling, with luminosities of ~ 3 — 10 x 10^^ erg/s, 
can remove the pressure support. Even for cases where the physical effects of hyperons were in- 
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eluded, whieh effeetively soften the EOS and reduce the maximum ahowed mass for an isolated NS 
to 1.8M0, the KT collaboration |259) still finds that thermal support can stabilize HMNS with 
masses up to 2.7 Mq. 

Using a Carpet/Cactus-based hydrodynamics code called Whisky [19 that works within the 
BSSN formalism (a version of which has been publicly released as GRHydro within the Einstein 
Toolkit [90]), the Whisky collaboration has analyzed the dependence of disk masses on binary 
parameters in some detail. For mass ratios q = 0.7 — 1.0 [236] . they found that bound disks with 
masses of up to 0.2 A/0 can be formed, with the disk mass following the approximate form 

A/disk = 0.039(Af„,ax - Aftot) + 1.115(1 - g)(AW - A/ft); A/„,ax = 1.139(1 + q)A'h , (34) 

where A/^nax the maximum mass of a binary system for a given EOS (F = 2 ideal gas for these 
calculations), A/^ is the maximum mass of an isolated non-rotating NS for the EOS, and A/tot the 
mass of the binary, with all masses here defined as baryonic. The evolution of the total rest mass 
present in the computational domain for a number of simulations is shown in Figure 1111 
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Figure 11: Evolution of the total rest mass Aftot of the remnant disk (outside the BH horizon) 
normalized to the initial value for NS-NS mergers using a F = 2 polytropic EOS with differing 
mass ratios and total masses (see Figure 5 of [236]). The order of magnitude of the mass fraction 
in the disk can be read off the logarithmic mass scale on the vertical axis. The curves referring to 
different models have been shifted in time to coincide at icon- 



6.4.2 Magnetized NS NS mergers 

Using the HAD code described in [6 that evolves the GHG system on an AMR-based grid with 
CENO reconstruction techniques, Anderson et al. [Q performed the first study of magnetic effects 
in full GR NS-NS mergers [7 . Beginning from spherical NSs with extremely strong poloidal 
magnetic fields (9.6 x lO^^G, as is found in magnetars), their merger simulations showed that 
magnetic repulsion can delay merger by 1-2 orbits and lead to the formation of magnetically 
buoyant cavities at the trailing end of each NS as contact is made (see Figure [T5)) . although 
the latter may be affected by the non-equilibrium initial data. Both effects would have been 
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greatly reduced if more realistic magnetic fields strengths had been considered. Magnetic fields 
in the HMNS remnant, which can be amplified through dynamo effects regardless of their initial 
strengths, helped to distribute angular momentum outward via the magneto-rotational instability 
(MRI), leading to a less differentially rotating velocity profile and a more axisymmetric remnant. 
The GW emission in the magnetized case was seen to occur at lower characteristic frequencies and 
amplitudes as a result. 

The UIUC group was among the first to produce fully self-consistent GRMHD results [86] . 
Using a newly developed Cactus-based code, they performed the first studies of unequal-mass 
magnetized NS-NS mergers [170] . Using poloidal, magnetar-level initial magnetic fields, Liu et al. 
found that magnetic effects are essentially negligible prior to merger, but can increase the mass 
in a disk around a newly formed BH moderately, from 1.3% to 1.8% of the total system mass 
for mass ratios of g = 0.85 and F = 2. They point out that MHD effects can efficiently channel 
outflows away from the system's center after collapse |290j , and may be important for the late-stage 
evolution of the system. 

In [115) . the Whisky group performed simulations of magnetized mergers with field strengths 
ranging from 10^^ to lO^^G. Agreeing with the UIUC work that magnetic field strengths would 
have essentially no effect on the GW emission during inspiral, they note that magnetic effects 
become significant for the HMNS, since differential rotation can amplify B-fields, with marked 
deviations in the GW spectrum appearing at frequencies of few ^ 2kIIz. They also point out that 
high-order MHD reconstruction schemes, such as third-order PPM, can produce significantly more 
accurate results that second-order limiter-based schemes. A follow-up paper [116j showed that a 
plausible way to detect the effect of physically realistic magnetic fields on the GW signal from a 
merger was through a significant shortening of the timescale for a HMNS to collapse, though a 
third-generation GW detector could perhaps observe differences in the kHz emission of the HMNS 
as well. 

More recently, they have used very long-term simulations to focus attention on the magnetic 
field strength and geometry found after the remnant collapses to a BH [237' They find that the 
large, turbulent magnetic fields {B ^ 10^^ G) present in the initial binary configuration are boosted 
exponentially in time up to a poloidal field of strength 10^^ G in the remnant disk, with the field 
lines maintaining a half-opening angle of 30° along the BH spin axis, a configuration thought to 
be extremely promising for producing a SGRB. The resulting evolution, shown in Figure I12[ is 
perhaps the most definitive result indicating that NS-NS mergers should produce SGRBs for some 
plausible range of initial parameters. 

It is worth noting that all magnetized NS-NS merger calculations that have been attempted to 
date have made use of unphysically large magnetic fields. This is not merely a convenience designed 
to enhance the role of magnetic effects during the merger, though it does have that effect. Rather, 
magnetic fields are boosted in HMNS remnants by the MRI, whose fastest growing unstable mode 
depends roughly linearly on the Alfven speed, and thus the magnetic field strength. In order to 
move to physically reasonable magnetic field values, one would have to resolve the HMNS at least 
a factor of 100 times better in each of three dimensions, which is beyond the capability of even the 
largest supercomputers at present, and likely will be for some time to come. 

6.4.3 GW emission 

In [144] , the KT collaboration found a nearly linear relationship between the GW spectrum cutoff 
frequency /cut and the NS compactness, independent of the EOS, as well as a relationship between 
the disk mass and the width of the kHz hump seen in the GW energy spectrum. While /cut is a 
somewhat crude measure of the NS compactness, it occurs at substantially lower frequencies than 
any emission process associated with merger remnants, and thus is the parameter most likely to 
be accessible to GW observations with a second generation detector. 
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Figure 12: Evolution of the density in a NS-NS merger, with magnetic field lines superposed, taken 
from Figure 1 of |237j . The first panel shows the binary shortly after contact, while the second 
shows the short-lived HMNS remnant shortly before it collapses. In the latter two panels, a BH 
has already formed, and the disk around it winds up the magnetic field to a poloidal geometry of 
extremely large strength, ^ lO^^G, with an half-opening angle of 30°, consistent with theoretical 
SGRB models. 

The qualitative form of the high-frequency components of the GW spectrum is primarily deter- 
mined by the type of remnant formed. In Figures and [T5l we show h{t) and h{f), respectively, 
for four of the runs calculated by the KT collaboration and described in [133] . Type I collapses 
are characterized by a rapid decrease in the GW amplitude immediately after the merger, yielding 
relatively low power at frequencies above the cutoff frequency. Type II and III mergers yield longer 
periods of GW emission after the merger, especially the latter, with the remnant oscillation modes 
leading to clear peaks at GW frequencies /gw = 2-4 khz that should someday be detectable by 
third generation detectors like the Einstein Telescope, or possibly even by advanced LIGO should 
the source be sufficiently close {D < 20Mpc) and the high-frequency peak of sufficiently high 
quality [260]. 

Using new multi-orbit simulations of NS-NS mergers, Baiotti et al. [TH [T5] showed that the 
semi-analytic effective one-body (EOB) formalism severely underestimates high-order relativistic 
corrections even when lowest-order finite-size tidal effects were included. As a result, phase errors of 
almost a quarter of a radian can develop, although these may be virtually eliminated by introducing 
a second-order "next-to-next-to-leading order" (NNLO) correction term and fixing the coefficient 
to match numerical results. The excellent agreement between pre-merger numerical waveforms and 
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Figure 13: Fluid density iscontours and magnetic field distribution (in a plane slightly above the 
equator) immediately after first contact for a magnetized merger simulation (see Figure 1 of [7|). 
The cavities at both trailing edges are attributed to magnetic pressure inducing buoyancy. 

the revised semi-analytic FOB approximant is shown in Figure 1161 
6.4.4 Binary eccentricity 

The effects of binary eccentricity on NS-NS mergers was recently studied by the Jena group [121] . 
Such systems, which would indicate dynamical formation processes rather than the long-term 
evolution of primordial binaries, evolve differently in several fundamental ways from binaries that 
merger from circular orbits. For nearly head-on collisions, they found prompt BH formation and 
negligible disk mass production, with only a single GW burst at frequencies comparable to the 
quasi-normal mode of the newly formed BH. For a collision in which mass transfer occurred at the 
first passage but two orbits were required to complete the merger and form a BH, a massive disk 
was formed, containing 8% of the total system mass, more than half of which was accreted by the 
BH within lOOM of its formation. Between the first close passage and the second, during which 
the two NS merged, the GW signal was seen to be quasi-periodic, and a a frequency comparable 
to the fundamental oscillation mode of the two NS, a result that was duplicated in a calculation 
for which the periastron fell outside the Roche limit and the eccentric binary survived for the full 
duration of the run, comprising several orbits. 
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Figure 14: Dimensionless GW strain Dh/mo, where D is the distance to the source and mo the total 
mass of the binary, versus time for four different NS-NS merger calculations, taken from Figures 5 
and 6 of [133) . The different merger types become apparent in the post- merger GW signal, clearly 
indicating how BH formation rapidly drives the GW signal down to negligible amplitudes. 



6.5 Simulations including microphysics 

In parallel to efforts in full GR, there has also been great progress in numerical simulations that 
include approximate relativistic treatments but a more detailed approach to microphysical issues. 
The first simulations to use a realistic EOS for NS-NS mergers were performed by Ruffert, Janka, 
and collaborators |2491I138|I250) . who assumed the Lattimer-Swesty EOS for their Newtonian PPM- 
based Eulerian calculations. They were able to determine a physically meaningful temperature for 
NS-NS merger remnants of 30-50 MeV, an overall neutrino luminosity of roughly lO^'^ erg/s for 
tens of milliseconds, and a corresponding annihilation rate of 2— 5x 10^° erg/sec given the computed 
annihilation efficiencies of a few parts in a thousand This resulted in an energy loss of 2—4 x 10^^ erg 
over the lifetime of the remnant [138) . a value later confirmed in multigrid simulations that replaced 
the newly formed HMNS by a Newtonian or quasi-relativistic BH surrounded by the bound material 
making up a disk |247j . The temperatures in the resulting neutron-rich {Yf. ^ 0.05 — 0.2) remnant 
were thought to be encouraging for the production of r- process elements [250] , although numerical 
resolution of the low-density ejecta limited the ability to make quantitatively accurate estimates of 
its exact chemical distribution. Further calculations, some of which involved unequal-mass binaries, 
indicated that the temperatures and electron fractions in the ejecta were likely not sufficient to 
produce solar abundances of r- process elements [248] . with electron fractions in particular smaller 
than those set by hand in the r-process production model that appears in [110II241J . More recently. 
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Figure 15: Effective strain at a distance of 100 Mpc shown as a function of the GW frequency 
(sohd red curve) for the same four merger calculations depicted in 1141 taken from Figures 5 and 
6 of |133j . Post-merger quasi-periodic oscillations are seen as broad peaks in the GW spectrum 
at frequencies /gw = 2-4 kHz. The blue curve shows the Taylor T4 result, which represents a 
particular method of deducing the signal from a 3PN evolution. The thick green dashed curve and 
orange dot-dashed curves depict the sensitivities of the second-generation Advanced LIGO and 
LCGT (Large Scale Cryogenic Gravitational Wave Telescope) detectors, respectively, while the 
maroon dashed curve shows the sensitivity of a hypothetical third-generation Einstein Telescope. 

it was suggested |122j that the decompression of matter originally located in the inner crust of a NS 
and ejected during a merger has a nearly solar elemental distribution for heavy r-process elements 
{A > 140). This indicates that NS-NS mergers may be the source of the observed cosmic r-process 
elements should there be sufficient mass loss per merger event , M^j ^ 3 — 5 x 10^^ Mq, although 
these amounts have yet to be observed in full GR simulations which have often admittedly been 
performed using cruder microphysical treatments. 

In [240] ■ Rosswog and Davies included a detailed neutrino leakage scheme in their calculations 
and also adopted the Shen EOS for several calculations, finding in a later paper [244] that the 
gamma-ray energy release is roughly IQ^^erg, in line with previous results from other groups, but 
noting that the values would be significantly higher if temperatures in the remnant were higher, 
since the neutrino luminosity scales like a very high power of the temperature. These calculations 
also identified NS-NS mergers as likely SGRB candidates given the favorable geometry j245) , and 
the possibility that the MRI in a HMNS remnant could dramatically boost magnetic fields on the 
sub-second timescales characterizing SGRBs |246] . Rosswog and Liebendorfer |242] found that 
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Figure 16: Comparison between numerical waveforms, shown as a solid black line, and semi- 
analytic NNLO EOB waveforms, shown as a red dashed line (top panel), taken from Figure 14 
of |15j . The top panels show the real parts of the EOB and numerical relativity waveforms, and 
the middle panels display the corresponding phase differences between waveforms generated with 
the two methods. There is excellent agreement between with the numerical waveform almost up 
to the time of the merger as shown by the match of the orbital frequencies (bottom panel) . 



electron antineutrinos v^. dominate the emission, as had Ruffert and Janka |248j , though the exact 
thermodynamic and nuclear profiles were found to be somewhat sensitive to the properties of 
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the EOS model. More recently, using the VULCAN 2-diniensional multi-group flux-limited-diffusion 
radiation hydrodynamics code |171j to evaluate slices taken from SPH calculations, Dessart et 
al. [52] found that neutrino heating of the remnant material will eject roughly 10"'* Mq from the 
system. 

Price and Rosswog |229| 1243] performed the first MHD simulation of merging NS-NS binaries 
using an SPH code that included magnetic field effects, finding that the Kelvin-Helmholtz unstable 
vortices formed at the contact surface between the two NSs could boost magnetic fields rapidly up 
to 10"G. These resuhs were not seen in GRMHD simulations, where gains in the magnetic field 
strength generated by dynamos were limited by the swamping of the vortex sheet at the surface 
of contact by rapidly infalling NS material that went on to form the eventual HMNS or BH [116] . 
Longer-term simulations did note that shearing instabilities were able to support power-law, or 
perhaps even exponential, growth of the magnetic fields on longer timescales 10 s of ms), which 
augurs well for NS-NS mergers as the central engines of SGRBs \2S7\ . 

An effort to identify potential observational differences between NSs and COs with quark- 
matter interiors has been led by Oechslin and collaborators. Using an SPH code with CF gravity, 
Oechslin et al. |206[ 1208] considered mergers of NSs with quark cores described by the MIT bag 
model j68 |ll02|ll41] . which have significantly smaller maximum masses than traditional NSs. They 
found the hybrid nuclear-quark EOS yielded higher ISCO frequencies for NSs with masses > 1.5 Mq 
and slightly larger GW oscillation frequencies for any resulting merger remnant compared to purely 
hadronic EOS. Bauswein et al. [34] followed up this work by investigating whether "strangelets" , 
or small lumps of strange quark matter, would be ejected in sufficient amounts throughout the 
interstellar medium to begin the phase transition that would convert traditional hadronic NSs into 
strange stars. They determined that the total rate of strange matter ejection in NS-NS mergers 
could be as much as 10~® Mq per year per galaxy or essentially zero depending on the parameters 
input into the MIT bag model, with the upper values clearly detectable by orbiting magnetic 
spectrometers such as the AMS-02 detector that was recently installed on the International Space 
Station [18H 1147] . Further calculations concluded that the mergers of strange stars produce a 
much more tenuous halo than traditional NS mergers, more rapid formation of a BH, and higher 
frequency ringdown emission [35] , as we show in Figure 1171 

Oechslin, Janka, and Marek also analyzed a wide range of EOS models using their CF SPH 
code, finding that matter in spiral arms was typically cold and that the dynamics of the disk formed 
around a post-merger BH depends on the initial temperature assumed for the pre-merger NS [205] . 
They also determined that the kHz GW emission peaks produced by HMNSs may help to constrain 
various parameters of the original NS EOS, especially it's high-density behavior ^204,, with further 
updates to the prediction provided by Bauswein and Janka [35]. Most recently, Stergioulas et 
al. [291] studied the effect of non-linear mode couplings in HMNS oscillations, leading to the 
prediction of a triplet peak of frequencies being present or low mass (Mns = 1-2 — 1.35 M©) 
systems in the kHz range. 

6.6 Comparison to BH— NS merger results 

While there is a history of Newtonian, quasi-relativistic, post-Newtonian, and CF gravitational 
formalisms being used to perform BH-NS merger simulations, their results are nearly always 
quantitatively, if not qualitatively, different than full GR simulations, and we focus here on the 
latter (see [279] for a more thorough historical review). Most of the groups that have performed full 
GR NS-NS merger calculations have also published results on BH-NS mergers, including Whisky 
(for head-on collisions) [l72], KT 1285, 284, 278,|327l|27ll[l53l[l52], UIUC [SHIMllM], HAD [67], 
as well as the SXS collaboration [BTj i84» il07t JOB] and Princeton (for elliptical mergers) [2891 155] . 
Summarizing the results of these works, we get a rather coherent picture, which we describe below. 
The GW signal from BH-NS mergers is somewhat "cleaner" than that from NS-NS mergers. 
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Figure 17: Evolution of a binary strange star merger performed using a CF SPH evolution, taken 
from Figure 4 of 35^ . The "spiral arms" representing mass loss through the outer Lagrange points 
of the system are substantially narrower than those typically seen in CF calculations of NS-NS 
mergers with typical nuclear EOS models. 



since the disruption of the NS and its accretion by the BH rapidly terminate the GW emission. 
In general, 3PN estimates model the signal well until tidal effects become important. The more 
compact the NS, the higher the dimensionless "cutoff frequency" Mtot/cut at which the GW energy 
spectrum plummets, with direct plunges in which the NS is swallowed whole typically yielding 
excess power near the frequency maximum from the final pre-merger burst. For increasingly 
prograde BH spins, there is more excess power over the 3PN prediction at lower frequencies, but 
also a lower cutoff frequency marking the plunge (see the discussion in [279] ). From an observational 
standpoint, the deviations from point-mass form become more visible for a higher mass BH- 
NS system, because frequencies scale characteristically like the inverse of the total mass. The 
distinction is particularly important for Advanced LIGO, as systems with Mbh ^ SA/q typically 
yield cutoff frequencies within the advanced LIGO band at source distances of 13 ~ 100 Mpc, while 
for lower-mass systems the cutoff occurs at or just above the upper end of the frequency band. 
This is significantly different than the situation for NS-NS mergers, in which the characteristic 
frequencies corresponding to the merger itself typically fall at frequencies above the advanced LIGO 
high-frequency sensitivity limit, and those corresponding to remnant oscillations in the range 2- 
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4kHz, which will prove a challenge even for third-generation GW detectors. 

Disk masses for BH-NS mergers were found to be extremely small in the first calculations, all 
performed using non-spinning BHs |285ll284ll91) . but have since been corrected to larger values once 
more sophisticated grid-based schemes and atmosphere treatments were added to those codes. More 
recent results indicate disk masses for reasonable physical parameters can be as large as 0.4 M©, for 
highly-spinning (oBH/Af = 0.9) mergers |107j . with values of 0.035 — 0.05 Mq characterizing non- 
spinning models with mass ratios g ss 1/5 [152j . Mass loss into a disk is suppressed by misaligned 
spins, especially for highly-inclined BHs, so the aligned cases should currently be interpreted as 
upper limits for the disk mass when alignment is varied |107) . Overall, disk masses for BH-NS 
merger remnants are comparable to those from NS-NS merger remnants, and may not be clearly 
distinguishable from them based solely on the emission properties of the disk. For BH-NS mergers 
with mass ratios g = 1/3 and prograde spins of dimensionless magnitude obh/A/ ~ 0.5, the disk 
parameters found after a run performed with the inclusion of a finite-temperature NS EOS |84j 
indicated that the neutrino luminosity from the disk might be as high as 10^"^ erg/s. While NS-NS 
merger simulations have led to predictions of neutrino luminosities a few times larger than this, the 
result does indicate that BH-NS mergers are also plausible SGRB progenitor candidates, possibly 
with lower characteristic luminosities than bursts resulting from NS-NS mergers. 

The role of magnetic fields in BH-NS mergers has only been investigated recently [STJ [M] , in 
simulations that apply an initially poloidal magnetic field to the NSs in the binary. Magnetic fields 
were found to have very little effect on the resulting GW signal and the mass accretion rate for the 
BH for physically reasonable magnetic field strengths, with visible divergences appearing only for 
B ^ lO^^G |94j, which is not particularly surprising. Just as in NS-NS mergers, magnetic fields 
play very little role during inspiral, and unlike the case of NS-NS mergers there is no opportunity 
to boost fields at a vortex sheet that forms when the binary makes contact, nor in a HMNS via 
differential rotation. While the MRI may be important in determining the thermal evolution and 
mass accretion rate in a post-merger disk, such effects will likely be observable primarily on longer 
timescales. 

Just as fuU-GR NS-NS simulations do not indicate that such mergers are likely sources of 
the r-process elements we observe in the universe, BH-NS simulations in full GR make the same 
prediction: no detectable mass loss from the system whatsoever, at least in the calculations per- 
formed to date. The picture may change when even larger prograde spins are modeled, since this 
should lead to maximal disk production, or if more detailed microphysical treatments indicate that 
a significant wind can be generated from either a HMNS or BH disk and unbind astrophysically 
interesting amounts of material, but neither has been seen in the numerical results to date. 

As is seen in NS-NS mergers, the pericenter distance plays a critical role in the evolution of 
eccentric BH-NS mergers as well. Large disk masses containing up to O.3M0, with an unbound 
fraction of roughly 0.15Mq, can occur when the periastron separation is located just outside the 
classical ISGO, with GW signals taking on the characteristic zoom- whirl form predicted for elliptical 
orbits [289'. In between pericenter passages, radial oscillations of the neutron star produce GW 
emission at frequencies corresponding to the f-mode for the NS as well [88]. 
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7 Summary and Likely Future Directions 



Returning to the questions posed in the previous section, we can now provide the current state of 
the field's best answers, though this remains a very active area of research and new results will 
certainly continue to modify this picture. 

1. With regard to the final fate of the merger remnant, calculations using full GR are required, 
but the details of the microphysics do not seem to play a very strong role. It is now possible 
to determine whether or not a pair of NSs with given parameters and specified EOS will form 
a BH or HMNS promptly after merger, and to estimate whether a HMNS will collapse on a 
dynamical timescale or one of the longer dissipative timescales (see, e.g., |133) ). For NS-NS 
binaries with sufficiently small masses, it is also possible to determine quickly whether the 
remnant mass is below the supramassive limit for which a NS is stabilized against collapse by 
uniform rotation alone, and thus would be unlikely to collapse, barring a significant amount of 
fallback accretion, unless pulsar emission or magnetic field coupling to the outer disk reduced 
the rotation rate below the critical value. This scenario likely applies only for mergers where 
the total system mass is relatively small Mtot ^ 2.5 — 2.6 M© [71], even taking into account 
the current maximum observed NS mass of Mns = 1.97 M© [ST]. Based on the wide arrays 
of EOS models already considered, it is entirely possible to infer the likely fate for sets of 
parameters and/or EOS models that have not yet been simulated, although no one has yet 
published a "master equation" that summarizes all of the current work into a single global 
form. While magnetic fields with realistic magnitudes are unlikely to affect the BH versus 
HMNS question |1701I116] . finite-temperature effects might play a nontrivial role should NSs 
be sufficiently hot prior to merger |260J (and see also [205] ). In the end, by the time the 
second generation of GW detectors make the first observations of mergers, the high-frequency 
shot-noise cutoff will prove to be a bigger obstacle to determining the fate of the remnant 
than any numerical uncertainty. A schematic diagram showing the possible final fates for 
a NS-NS merger along with the potential EM emission (see Figure 21 of [277] ) is shown in 
Figure m 

2. GW emission during merger is also well-understood, though there are a few gaps that need to 
be filled, with full GR again a vital requirement. While the FN inspiral signal prior to merger 
is very well understood, finite-size tidal effects introduce complications beyond those seen in 
BH-BH mergers, yet the longest calculations performed to date [15] encompass fewer orbits 
prior to merger than the longest BH-BH runs |257j . As noted in [15] and elsewhere, longer 
calculations will likely appear over time, helping to refine the prediction for the NS-NS merger 
GW signal as the binary transitions from a FN phase into one that can only be simulated 
using full GR, and teasing out the NS physics encoded in the GW signal. It seems clear from 
the published work that the emission during the onset of the merger is well-understood, as 
is the very rapid decay that occurs once the remnant collapses to a BH, either promptly or 
following some delay. GW emission from HMNSs has been investigated widely, and there have 
been correlations established between properties of the initial binary and the late-stage high- 
frequency emission (see, e.g., [1441 1116j ). but given that magnetic fields, neutrino cooling, 
and other microphysical effects seem to be important, a great deal of work remains to be 
done. Ferhaps more importantly, since HMNSs emit radiation at frequencies well beyond 
the shot-noise limit of even second-generation GW detectors, while the final inspiral occurs 
near peak sensitivity, it is likely that the first observations of NSs will constrain the nuclear 
EOS (or perhaps the quark matter EOS [35] ) primarily via the detection of small finite-size 
effects during inspiral. Since QE calculations are computationally inexpensive compared to 
numerical merger simulations, there should be much more numerical data available about the 
inspiral stage than other phases of NS-NS mergers, which should help optimize the inferences 
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Figure 18: Summary of potential outcomes from NS-NS mergers, taken from |277| . Here Mthr is 
the threshold mass (given the EOS) for collapse of a HMNS to a BH, and Qm is the binary mass 
ratio. 'Small', 'massive', and 'heavy' disks imply total disk masses Mdisk ^ 0.01 M©, 0.01 < 
A^disk ^ 0.03 Mq, and A/disk ^ 0.05 M0, respectively. 'B-field' and 'J-transport' indicate potential 
mechanisms for the HMNS to eventually lose its differential rotation support and collapse: magnetic 
damping and angular momentum transport outward into the disk. Spheroids are likely formed only 
for the APR and other stiff EOS models that can support remnants with relatively low rotational 
kinetic energies against collapse. 



to be drawn from future observations. 

3. Determining the mass of the thick disk that forms around a NS-NS merger remnant remains a 
very difficult challenge, since its density is much lower and harder to resolve using either grid- 
based or particle-based simulations. The parameterization given by Eq. [M] [236] is generally 
consistent with the GR calculations of other groups (see, e.g., |133] ). and seems to reflect 
a current consensus. It is also clear that disk masses around HMNSs (up to 0.2 A/q) are 
significantly larger than those forming around prompt collapses, which are limited to about 
0.05 Mq. 

4. It is likely that several orders of magnitude more mass-energy are present in the remnant 
disk than is observed in EM radiation from a SGRB. Modeling the emission from the disk 
(and possibly a HMNS) remains extremely challenging. Neutrino leakage schemes have been 
applied in both approximate relativistic calculations |2481 1242] and full GR |260J , and a more 
complex flux-limited diffusion scheme has been applied to the former as a post-processing 
step [82] , but there are no calculations that follow in detail the neutrinos as they flow outward, 
annihilate, and produce observable EM emission. At present, nuclear reactions are typically 
not followed in detail; rather, the electron fraction of the nuclear material, Y^, is evolved, 
and used to calculate neutrino emission and absorption rates. 

5. Magnetic fields, on the other hand, are starting to be much better understood. B-fields do 
seem to grow quite large through winding effects, even during the limited amount of physical 
time that can currently be modeled numerically [71 1170| 13271 1237] , with some calculations 
indicating exponential growth rates. The resulting geometries seem likely to produce the 
disk/jet structure observed throughout astrophysics when magnetized objects accrete mate- 
rial, which span scales from stellar BHs or pre-main sequence stars all the way up to active 
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galactic nuclei |237j . 



6. While recent numerical simulations have strengthened the case for NS-NS mergers as SGRB 
progenitors, full GR calculations have not generated much support for the same events yield- 
ing significant amounts of r-process elements. Noting the standard caveat that low-density 
ejecta are difficult to model, there is still tension between CF calculations producing ejecta 
with the proper temperatures and masses to reproduce the observed cosmic r-process abun- 
dances (see, e.g., |122J ). and full GR calculations that produce almost no measurable unbound 
material whatsoever. 

While NS-NS merger calculations have seen tremendous progress in the past decade, the future 
remains extremely exciting. Between the addition of more accurate and realistic physical treat- 
ments, the exploration of the full phase space of models, and the linking of numerical relativity to 
astrophysical observations and GW detection, there remain many unsolved problems that will be 
attacked over the course of the next decade and beyond. 
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A Field evolution equations 

In the BSSN evolution system, we define the following variables in terms of the standard ADM 
4-metric gij, 3- metric 7,j, and extrinsic curvature Kij: 



log [det7y] = log tjj, 



K 
pi 
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The evolution system consists of 15 equations for the various field terms, 
doK = 



do(t> 
dolij 
doAij 



-40 



aRij + ai?f • - DiD^a 



TF 



2 ~. 



+ aKAij - 2aAikA) + 2Ak(idj)p^ - -1^,9^/3'= - Snae-^^Sf/ 



-2A'^dia + 2a 



t\iA'''+6A^^djcp--rKj 



(35) 

(36) 
(37) 

(38) 

(39) 

(40) 
(41) 

(42) 

(43) 
(44) 



where the matter source terms contain various projections of the stress-energy tensor, defined 
through the relations 



E 



Si 



1 



n' 



"n'^T^, = ^ (Too - 2/3'Toi + p'p^Tij) , 
-n^^j-T^, = -^{To^-(3^Tij), 



Si' 



(45) 

(46) 
(47) 



Wc have introduced the notation da = di — [3^ dj . All quantities with a tilde involve the conformal 
3-metric 7^ , which is used to raise and lower indices. In particular, Di and r*j refer to the covariant 
derivative and the Christoffcil symbols with respect to 7^^. Parentheses indicate symmetrization 
of indices, and the expression [■■■]'^^ denotes the trace- free part of the expression inside the 
brackets. In the BSSN approach, the Ricci tensor is typically split into two pieces, whose respective 
contributions are given by 



Ri 



Rl 



-2D,D,(, 



2%D''Dk^ + AD,^ Dj(f, ~ A%D''(I) Dk^ 



(48) 
(49) 



These equations must be supplemented with gauge conditions that determine the evolution of 
the lapse function a and shift vector Noting that some groups introduce slight variants of these, 
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the moving puncture gauge conditions that have become popular for all GR merger calculations 
involving BHs and NSs typically take the form 

doa = -2aK, (50) 

dtP' = (51) 

dtB' = dtr-7-iB\ (52) 

where is an intermediate quantity used to convert the second-order "Gamma-driver" shift 
condition into a pair of first-order equations, and 77 is a user-prescribed term used to control 
dissipation in the simulation. Note that it is possible to replace the three instances of dt in 
the shift evolution equations Eqs. [5T] and [52] bv the shift-advected time derivative {dt — P-'dj) 
without changing the stability or hyperbolicity properties of the evolution scheme; in both cases 
moving punctures translate smoothly across a grid over long time periods |317| and both systems 
are strongly hyperbolic so long as the shift vector does not grow too large within the simulation 
domain [128] . 

The generalized harmonic formulation involves recasting the Einstein field equations, Eq. [5]in 
the form 

Rf^y = Stt - ]^g^,vT^ (53) 
and, after some tensor algebra, rewriting the Ricci tensor in the form 

-2i?^, = g^'g^p,^s + g-'^pgcsn + g^dps^-i + 2i^(a,« - 2i/,-ri^ + 2r]^r^„ (54) 
The Christoffcl coefficients are calculated from the full 4-metric, 

1 

2' 



r^7 = ^3"'' [gpsn + 9i&:fi - gpiA (55) 



and the gauge source terms i?^ are defined in Eq. [25] 

Given well-posed initial data for the metric and its first time derivative (since the system is 
second-order in time according to Eq. [5^ . the evolution of the system may be treated by a first- 
order reduction that specifies the evolution of the four functions if along with the spacetime 
metric g^j/, its projected time derivatives Ti^n = —n"dag^u, and its spatial derivatives 

'^ifiy = d.^g^,^ (56) 

subject to a constraint specifying that the derivative terms remain consistent with the metric 
g^i^ in time. In practice, one typically introduces a constraint for the source functions, defining 

C^" = H^' - Dx^ (57) 

and then modifies the evolution equation by appending a constraint damping term to the RHS of 
the stress energy-term (following [1271 1227] 

-87r(2r^^ - g^„T) - 8tt{2T^^ - g^^^T) - K{n^C^ + n^C^ - g^^^n^Cs) (58) 

where is the unit normal vector to the hypersurface (see Eq. I15p . The gauge conditions used in 
the first successful simulations of merging BH binaries |227] consisted of the set 

OHt = -6^^ + = (59) 

with £,i,^2,V constant parameters used to tune the evolution. The first one drives the coordinates 
toward the ADM form and the latter provides dissipation. The binary NS-NS merger work of [7] 
chose harmonic coordinates with = 0. 
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B GR Hydrodynamical and MHD equations 

In what follows, we will adopt the stress energy tensor of an ideal relativistic fluid, 

T^"^ = phu^'u'' + g^'^P , (60) 

where p, P, and m'' are the rest mass density, pressure, and fluid 4-velocity, respectively, and 

h = l + e + P/p (61) 

is the relativistic specific enthalpy, with e the specific internal energy of the fiuid. 

The equations of ideal GR hydrodynamics [184] may be derived from the local GR conservation 
laws of mass and energy-momentum: 



V. = 0, V.T^'' = , 



(62) 



where denotes the covariant derivative with respect to the 4- metric, and = pu^ is the mass 
current. 

The 3- velocity -y' can be calculated in the form 



where 



W = au° = {l- v'v,)-^/^ 
is the Lorentz factor. The contravariant 4-velocity is then given by: 



W 



— , u' = W v'-^ 



and the covariant 4-velocity is: 



uo = W{v'P, - a) 



Ui = Wvi 



(63) 
(64) 

(65) 
(66) 



To cast the equations of GR hydrodynamics as a first-order hyperbolic flux-conservative system 
for the conserved variables D, S^, and r, defined in terms of the primitive variables p, e, w% we define 



S' = y/^phW^v\ 
T = ^{phW^-P)-D, 

where 7 is the determinant of jij . The evolution system then becomes 

dV dF' „ 



with 



dt dx^ 

U = [D,S,,t], 



(67) 
(68) 
(69) 

(70) 



(71) 
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Here, v'^ = — /a. and Vj^^ are the 4-Christoffel symbols. 

Magnetic fields may be included in the formalism, in the ideal MHD limit under which we 
assume infinite conductivity, by adding three new evolution equations and modifying those above 
to include magnetic stress-energy contributions of the form 



(72) 



where the magnetic field seen by a comoving observer, ft** is defined in terms of the dual faraday 
tensor *f by the condition 

b'^=*F"i'u^ (73) 

where 6^ = V^b^ represents twice the magnetic pressure. With magnetic terms included, we may 
rewrite the stress-energy tensor in a familiar form by introducing magnetically modified pressure 
and enthalpy contributions: 



T^"" = ph*u>'u'' + P*g>"'; h* = l + e+?^^, P* 



9 

^+2 



(74) 

(75) 
(76) 

(77) 

and the conserved variable = y^S*, which are related to the comoving magnetic field 4- vector 
6'* through the relations 

WB'^Vk 



and redefine the conserved momentum and energy variables S'^ and r accordingly: 

= ^ph*W'^v' - ab%' 
T = ^{ph*W^-P*-{ab°f)-D 
Defining the (primitive) magnetic field 3- vector as 



a r 



*T7lOi 



6° 

b' 

bi 
b^ 



a 



B 



W 

B. 
W 

B^Bj 



W{B''vk)vi 
+ [B'vif, 



(78) 

(79) 

(80) 

(81) 
(82) 



we may rewrite the conservative evolution scheme in the form 



U = 



a X 



S = 



a X 



SjV' + SiP* -bjB'/W, 
TV' + P*v' - ab^B'/W, 

0, 




(83) 



(84) 



where the magnetic field evolution equation is just the rclativistic version of the induction equation. 
An external mechanism to enforce the divergence-free nature of the magnetic field, diB' must also 
be applied. 
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